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INTRODUCTION 

The GSFC Faculty Fellowship Program is a cooperative effort between Goddard, the 
American Society of Engineering Education (ASEE), the University of Maryland, and the 
Catholic University of America. 

Under this program university faculty members come to Goddard for a period often 
weeks during which they engage actively in research and at the same time participate in 
seminars related to their research. 

The objectives of this program are many; some of them are listed below: 

1. Stimulation of schools to become interested in the research problems confront- 
ing Goddard. 

2. Creation of interest on the part of the university faculty to continue their research 
after completing the formal program. 

3. Stimulation of our people professionally through associations with the university 
faculty and through participation in the program's seminars. 

4. Establishment of closer ties with the Universities. 

This report marks the completion of the sixth year of this program at Goddard, and 
represents the progress in one of its objectives. 



Director 
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PRE-VIBRATION MAGNETIC TESTING OF THE SAS-A 


Dr. Virgilio Acosta 
Associate Professor of Physics 
United States Naval Academy 
Annapolis, Maryland 21402 

1. BRIEF DESCRIPTION OF THE MAGNETIC TEST SITE FACILITY 

The Magnetic Test Site is a complex of eleven different buildings outside the main 
campus of GSFC and located in a wooded area near Good Luck Road, in Greenbelt, 
Maryland. A sketch of the distribution of the buildings, not on scale, appears in Figure 1. 



MAIN OBJECTIVES - 

The Magnetic Test Site consists of two principal test facilities; the Spacecraft Mag- 
netic Testing Facility (SMTF) and the Magnetic Instrument Test Lab (MITL). 

The SMTF is intended to furnish the GSFC with enough capability for the evaluation 
of some aspects relative to the control of the orientation of spacecrafts in a magnetic 
environment, similar to the one they are going to find in outer space. More specifically, 
the SMTF is intended to provide a controllable large magnetic field (maximum value 60K 
gamma) and spacewise large enough to accommodate spacecrafts and components bigger 
than those that can be tested in the MITL. Other capacilities include: 
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a. Controllable ways for the simulation of the relative motion of the spacecraft under 
test and the magnetic field in the outerspace in order to find the magnetic forces and 
torques acting on a flyii^ spacecraft. 

b. Mapping contours of the spacecraft’s intrinsic magnetic field in a controllable ex- 
ternal magnetic field. 

c. Study the effects of external magnetic forces and torques on the different subsys- 
tems and components of the spacecraft with the main purpose of having attitude con- 
trol or remove interference with magnetic instruments outboard the spacecraft. The 
SMTF is located in Building 4. 


DESCRIPTION OF FACILITIES 

Buildii^ 2 - The central office is located in this building as well as a machine shop 
and the controls of air conditioning in the other buildings. 

Building 3 - Contains a tri-axial Helmholtz's coil assembly consisting of 12 coils, 4 
horizontal, 4 oriented with the common axis in the NS direction, and 4 with the com- 
mon axis oriented in the EW direction. Maximum diameter of the coils is 22 ft. 
Magnetic testing of small spacecrafts or of sounding rockets can be done in this 
building. Maximum field in any direction is 60K gamma ± 0.2 gamma. The mag- 
netometer that is used in this facility is a triaxial flux-gate, maximum sensitivity 
is 0.1 gamma. Field can be rotated in any direction at a maximum rate of 100 
rad/s. The building is also equipped with a Rubidium vapor magnetometer. Maxi- 
mum sensitivity is 0.1 gamma. 

Building 5 - In this building a three-axis Braunbeck coil system (Figure 2) consist- 
ing of 12 coils is installed as in Buildii^ 3; 4 coils are horizontal, 4 with the common 
axis oriented in the NS direction and 4 in the EW direction. The main characteristics 
are maximum diameter 42 feet; completely non-magnetic structure; magnetic field 
0 - ± 60K gamma in any direction; accuracy 0.20 gamma; rotating field 0 - 60K gamma 
in any direction at a maximum rate of 100 rad/s. Large spacecrafts or rockets are 
tested in this building. The spacecraft magnetic test facilities are located in this 
building. 

Building 4 - The electronics and instrumentation for the coils in Buildings 3 and 5 
are located here. For data acquisition there are four fluxgate triaxial magnetometers 
with a maximum sensitivity of 0.1 gamma. Several three-axes flux-gate magnetom- 
eters (maximum sensitivity 0.1 gamma) and a Rubidium vapor magnetometer (maxi- 
mum sensitivity 0.1 gamma) are part of the associated equipment. 

Building 1 - The main facility in this building is a square coil assembly, 14 feet on 
each side. The main purpose is to produce a zero field at the center of the facility 
in order to cancel out the earth's magnetic field within 0.5 gamma in any direction. 

The facility contains the necessary instrumentation (three axes flux-gate magnetom- 
eter, voltmeter, printer, etc.) for performing magnetic field measurements. Com- 
ponents, subsystem testing and dipole determinations can be made in this facility. 

Building 8 - Magnetic Shelter. This contains the proper control for canceling out 
the diurnal fluctuations of the earth magnetic field within 0.5 gamma at the center 
of the coils in Building 3. 

Building 9 - Here are located the necessary electronic hardware to operate the 
controls in Building 8. 
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SPACECRAFT MAGNETIC 



Figure 2. 40 Foot Brounbeck Coil Assembly 




Building 10 - Magnetic Shelter. Contains the proper controls for cancelling out the 
diurnal fluctuations of the earth magnetic field within 0.5 gamma in a spherical 
volume about 6 feet in diameter with its center coinciding with the center of the 
Braunbeck coil system in Building 5. 

Building 11 - Here are located the necessary electronic hardware to operate the 
controls in Building 10. 

Building 6 - Quiet Lab 1. This facility is operated mainly by the personnel of the 
Space Science Directorate. The main facility in this building is a Braunbeck 12 foot 
coil assembly, consisting of 10 coils; 4 are horizontal, 4 are oriented with the com- 
mon axis in the NS direction, and 2 are oriented with the common axis pointing in 
the EW direction. It is possible to produce a zero field at the center of the coils 
and cancel out the diurnal fluctuations of the earth's magnetic field within 0.5 gamma 
in each direction. The building is equipped with several three axes flxix-gate mag- 
netometers (maximum sensitivity 0.1 gamma). The main purpose of the facility is 
magnetometer development. 

Building 7 - Quiet Lab 2. This facility, as the one installed in Building 6, is operated 
mainly by personnel of the Space Science Directorate. The main facility in this 
building is a cubical coil assembly, 15 feet on each side. It is possible to cancel 
out the earth's magnetic field at the center of the coils (within 0.5 gamma) in any 
direction and at the same time cancel out the diurnal fluctuations of the geomagnetic 
field (within 0.5 gamma). The building is equipped with several three axes flux-gate 
magnetometers with a maximum sensitivity of 0.1 gamma. The facility is mainly 
used in magnetometer development. 


2. DESCRIPTION OF THE PREPARATORY WORK DONE (INCLUDING DATA) 

FOR THE SAS-A PRE -VIBRATION MAGNETIC TEST. 

As one of the main goals of the actual testing of the spacecraft, the determination of 
its associated dipole- moment, it is necessary to simulate the geomagnetic field that the 
satellite is going to find when it is in flight in order to compute the magnetic torque that 
will be actir^ on it in order to have ground control on its attitude (orientation) and spin 
rate. The magnetic torque is computed by the equation: 


L =Mx B 


( 1 ) 


in which 

M is the magnetic moment in pole- cm. 

B is the magnetic field in oersted. 

L is the magnetic torque in dyne-cm. 

In order to measure the torque a very ii^enious device called a torquemeter was 
designed at the Magnetic Test Site. (Figure 3). The torquemeter that is in use at present 
is an improved version of the one represented in Figure 3 and was designed by Joseph 
C. Boyle of the staff of the Magnetic Test Site and in charge of the experimental testing 
done on spacecrafts, components and subsystems. 
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Figure 3. Schematics of the Torquemeter 


This improved torquemeter is capable of measuring torques in the vertical direction 
in the range 0.3 dyne- cm to 10® dyne-cm and can be loaded with a maximum load of 800 
pounds. In order to produce a known magnetic moment for calibrating the torquemeter 
and for nulling spacecraft, torques, two air core solenoids at right angles and in a hori- 
zontal plane are attached to the torquemeter. These air core coil solenoids can be ener- 
gized by means of a D.C. current and they have independent circuits, so they can be op- 
erated either separately or simultaneously. The purpose of these solenoids is to associ- 
ate to the torquemeter a magnetic dipole moment that can be varied at will because the 
other components of the meter are made of non- magnetic materials. Now by placing the 
torquemeter in an environmental magnetic field, a magnetic torque will be acting on it, 
the vertical components of which can be measured by remote control. As the two air- 
core solenoids are mounted at right angles and in a horizontal plane it is necessary to 
calibrate them in order to know the magnetic moment in the horizontal plane associated 
with the torquemeter. The calibration constant in pole-cm/amp was measured mag- 
netically by use of the torquemeter. In order to distinguish the air-core solenoids, we 
will designate them as the radial and tangential coils. If M is the magnetic moment as- 
sociated with one solenoid and I the energizing current, the calibration constant will be 
obtained by the expression; 


K= — ploe-cm/ amp. 


( 2 ) 


Magnetic Calibration of the Air-Coils of the Torquemeter - This was done in Build- 
ii^ 1, where one of the cubical coil systems is located (See Figure 3 for the experimental 
setup). 



Figure 4. Schematic Set-Up for the Static 
Air-Coil Calibration 


A represents a coil placed horizontally at the center of the cubical coil assembly 
with its axis pointing in the NS direction and B aligned in the same direction is the probe 
of a flux-gate magnetometer, the distance r was set at 6 feet. A zero magnetic field was 
set at the center of the coils, a regulated D.C. current of 1 amp was sent through the 
coil at A and the corresponding radial magnetic field B^ was measured with the magnetometer. 
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Then the magnetic moment M associated with the coil was computed using the equation: 



( 3 ) 


Several measurements were made with each coil and the average taken for M. The re- 
sults were as follows: 


Radial Coil K 


M R 
I 


2349 


Pole - cm 
Amp . 


Tangen tial 


Coil 


M t 
I 


1958 


Pole - cm 
Amp. 


Before proceeding to check dynamically the calibration of the two air-core coils of the 
torquemeter, we performed several preliminary tests relative to the torquemeter itself. 
In order to conduct these tests we moved the torquemeter to Building 5 where the 
Braunbeck coil system is installed, and three balance- noise tests were conducted, with 
the torquemeter unloaded and then loaded with about 100 pounds and finally with 400 
pounds. The reason for loading the torquemeter with 400 pounds was to simulate the 
approximate weight and moment of inertia of the SAS-A spacecraft that we were going to 
test later on. In this test the torquemeter was placed at the center of the Braunbeck coil 
system as it is shown in Figure 5 with the Z axis in a vertical position, D, and Dj are 
two flat squeeze film dampers and by pouring water in both they exert a damping effect 
on the vibrations of the torquemeter around the Z axis. The main goals of these balance 
noise tests were to bring the center of mass of the torquemeter to same point along the 
Z axis and to minimize its output noise by balancing the torquemeter once it is leveled. 



At this point it is convenient to mention that the output mechanical signal of the torquemeter 
(the torque acting on it) is transformed into an electrical signal by means of an air differ- 
ential capacitor (attached to the torquemeter) whose capacitance changes with the applied 
torque and torque fluctuations acting on it. The electrical output of the torquemeter is 
transmitted to a two channel recorder where it can be measvired and recorded. A block 
diagram of the setup used for the balance- noise is shown in Figure 6. 
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Torquemeter Filter Recorder 

Figure 6. Block Diagram of the Electric System 
used for the Balance-Noise Tests 


The object of the filter is to get rid of undesirable frequencies in the recorder input. 
The filter has three modes of operations: low pass, band pass, and high pass. Noise 
being recorded by the recorder come from different sources: natural seismic vibration 
of the building, vibrations produced by flyii^ air planes or any other moving vehicle close 
to the building, external or internal air currents, thunderstorms, etc. In order to mini- 
mize air currents inside the building, the Braunbeck coil system is provided with a poly- 
thylene covering (shroud) that can be placed around the instrument that is under test. As 
we said before, the main goal of this test was to reduce to a minimum the noise coming 
from the torque meter due to the above mentioned factors, but mainly to the lack of bal- 
ance of the torquemeter. As a result of different tests the noise was reduced to a mini- 
mum by adding two weights C and D of about 5 pounds each, attached to the upper moving 
platform of the torquemeter (Figure 5), both dampers full with water and at a band pass 
frequency of 1.00 Hz as input to the recorder. The test was conducted with no current in 
the Braunbeck assembly and without energizii^ the radial and tangential coils and with an 
applied torque of 681 dyne-cm applied to the torquemeter. This calibrating torque can be 
applied to the torquemeter by a very delicate balance (bell- crank) in which calibrated 
weights produced known torques acting along the Z axis of the torquemeter. 


SIGNAL TO NOISE RATIO 

As it is known, the vibrations of a body under the action of a driving sinusoidal torque 
T = Tg sin wt obey the differential equation 


Ix + cx + Kx = Tp sincot 


( 3 ) 


which is called the differential equation of a forced rotational oscillation. In this equation 
I = Moment of inertia around the axis of rotation 
C = Damping Factor 
K = Torsional constant 
X = Angular deflexion 

’’’ = in wt = applied sinusoidal torque along the axis of rotation. 

The corresponding ar^lar deflections are given by 



( 4 ) 
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in which 


Xst = constant static deflection under a constant torque and given by Hooke's Law 
T = in which t is the applied torque and K the torsional constant. 

f = Forcing frequency 

fg = Natural frequency 

c = Actual damping factor 

Cg = Critical damping factor 

The ratio: 



( 5 ) 


is called the magnification factor, a family of frequency response curves in which MF is 
plotted versus f/fg and c/Cg is taken as parameter, is shown in Figure 7. 

In order to find the signal to noise ratio, the following procedures was followed: 

a. Torquemeter balanced and loaded with 400 pounds, 

b. Dampers full 

c. Calibrating mechanical weight producing a Z axis torque of 681 dyne-cm, 

d. A 2 amp current was set in the EW (radial coil) in order to produce a magnetic 
dipole moment associated with the torquemeter and located in a horizontal plane. 


M.F. 



Figure 7. Frequency response curves for different values of c/cq. 
a) shows no damping c/c^ = 0, d) shows high damping, b) corre- 
sponds to the maximum M.F. at resonance frequency. 


9 




e. No magnetic field sets up by the Braunbeck assembly. 

f. Geomagnetic field of the earth present. 

Under these conditions the current in the E-W coil was suddenly reduced to zero and 
the output signal was recorded in the dual recorder. A typical curve obtained from the 
recorder appears in Figure 8. 

The main purpose of this operation was to obtain a decay curve from which we ob- 
tained the natural frequency as well as the damping c/Cq . This allowed us to "fictionalize" 
the response to an imaginary input torque signal. This fictional response was then com- 
pared to the actual measured noise with the filter set in the band pass mode at the fre- 
quencies tabulated (Table I) to obtain the signal to noise ratio. 

Several curves were obtained for different frequency inputs to the recorder as set by 
the filter. Results appear in Table I. 


Table I 


Torquemeter Responses for Different Frequencies Assuming Torquemeter 
signal = 1 for f = o and setting arbitrarily X^t =1 


f(Hz) 

Signal X 
with X St = 1 

Peak 

Noise from Recorder 

Signal 

Noise 

Relative Signal 
Noise 

0.25 

4,20 

15 

0.28 


0.283 (fo) 

10.00 

17 

0.59 

2.0 

0.50 

0.47 

8 

0.059 

0,2 

0.75 

0.17 

3 

0.057 

0.2 

1.00 

0.09 

2 

0.045 

0.16 

1.25 

0.05 

3 

0.017 

0.06 

1.50 

0.04 

3 

0.013 

0.05 


As can be seen, the maximum M.F. = 10 and corresponds to the natural frequency 
0.283 Hz of the loaded torquemeter. The natural frequency of oscillation of the torque- 
meter loaded with 400 pounds was determined experiment^ly by an analysis of the data 
obtained from the recorder. The values X/Xst were computed by using the equation: 
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(5) 




In order to find the relative damping c/c,, (which is necessary to compute the values 
of X/X^f according to Eg 5, the ratio c/cq must be calculated first. This was done in the 
followii^ way. It can be proved that if X_^ and X^+i (see Figure 8) are two successive 
peaks in the torquemeter output the following relationship follows for small damping: 


hence: 



C 



1 

2 7T 



(6) 


Several values of c/c^ were computed yielding an average: 


— = 0.05 

c„ 


The M.F. at resonance is then given by the equation: 



(7) 


Replacing the value c/cg = 0.05 in equation 7 the result coincides with the result of Table 
I (second row of table). 


Mechanical Calibration of Torquemeter and Correlation with both Air-Coil Dipole 
moments. Dynamic Calibration of the Air-Coils and Checking their Alignment. 

These tests were done in several steps as described here. The first test was done 
under the following conditions (Figure 5): 

a. Torquemeter setup at center of the Braunbeck coil assembly with the Z axis in a 
vertical position and the radial coil oriented east-west. 

b. No current in the Braunbeck coil assembly. 

c. Geomagnetics field of the earth present. 

d. Applied torques to the torquemeter by means of bell crank device: 681 dyne-cm. 

e. Filter in the low pass mode of operation at a cutoff frequency of 1 Hz. 

f. Energizing current in the air-coils attached to the torquemeter was 2 amp. 

g. Torquemeter load: 400 pounds. 
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The block diagram of the electrical setup used for this calibration appears in Figure 6. 
The calibration constant as recorded (average calibration) by the recorder was 32.7 
(dyne-cm/div.) at 1 volt/cm sensitivity in the recorder. This calibration was a me- 
chanical calibration of the torquemeter. 

The following step was a check of the proper alignment of the air-coils attached to 
the torquemeter to within 1 degree. 

Test was conducted with the torquemeter set as shown in Figure 5. Conditions for 
this test were as follows (Figure 5): 

a. Torquemeter set up at center of Braunbeck coil system and with the Z axis in a 
vertical position. 

b. Dampers full 

c. Torquemeter load: 400 pounds. 

d. No energizing current in the air-coils attached to the torquemeter. 

e. Electrical setup as in Figure 6. 

The alignment of the radial coil was performed first under a fixed field of 60K 
gamma, E-W direction provided by Braunbeck coil assembly. The result was acceptable. 
The criterion for coil alignment was that the torque produced when the radial coil was 
energized in the presence of an E-W field does not exceed that due to 1 degree misalign- 
ment. The maximum allowance for a 1 degree misalignment of the radial coil was 40 
pole-cm. 

At this time we also verified the coil constants of both coils (within 5% of the mag- 
netically obtained values) by measuring the torque produced in a known field applied at 
right angles with a known coil current. 

The same alignment test was then performed for the tangential coil under the same 
conditions as before but with a field of 60K gamma in the N-S direction. The result of 
this test was also considered as acceptable. 

The next step was a dynamic calibration of the torquemeter and correlation with the 
current in the air- core coils and at the same time a dynamic check on the proper align- 
ment of the air coils. 

The setup for the torquemeter in this step is again the same one shown in Figure 5. 

The electrical setup is shown in Figure 6. 

Test was conducted under the following conditions: 

a. Torquemeter setup at center of Braunbeck coil system and with Z axis in a 
vertical position. 

b. Dampers full 

c. Torquemeter load: 400 pounds. 

d. Current in the air- coils: 2 amp. 

e. Filter in the band-pass mode at 0.286 Hz (natural frequency). 
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The first coil tested was the tangential coil (cal. const. 1958 pole- cm/ amp). 

An oscillating field (maximum amplitude 60K gamma) oriented in the N-S direction 
and with a frequency of oscillation of 0.286 Hz was set up by the Braunbeck coil assembly. 

The result of this alignment check was: 41 dyne-cm/Amp (approximately 2° mis- 
alignment) for the tangential coil. The radial coil was not tested. 

The next step was a torque calibration of the air- coils attached to the torquemeter 
to check the magnetic calibration performed before. 

The results of the magnetic calibration were: 

Radial Coll = 2349 pole/cm/ amp 

Tangential Coil Kt = 1958 pole/cm/ amp 

In order to perform these torquemeter calibrations, the following conditions were 
followed: 

a. Torquemeter setup as in Figure 5. 

b. Torquemeter load: 400 pounds. 

c. Dampers full 

d. Electrical setup as in Figure 6. 

The tangential coil was tested first and a 60K gamma (E-W) field was set using the 
Braunbeck coil assembly. Then the tangential coil was energized with a current of 1 amp 
and the corresponding deflection in the recorder was measured. 

Same procedure was followed in the testii^ of the radial coil, but with a field of 30K 
gamma in the N-S direction. 

Results from these calibrations were as follows: 

a. Radial Coil K, = 2380 pole/cm/amp 

b. Tangential Coil = 1900 pole-cm/amp 

The results were considered as acceptable when compared with the former values 
obtained followir^ the magnetic procedure. 

Now a dynamic test was performed whose main purpose was to confirm that the 
best operational frequency for the SAS-A test was the natural frequency (0.283 Hz). 

A block diagram of the electrical setup is snown in Figure 9. 

The experimental setup is shown in Figure 5 but adding a fluxgate magnetometer near 
the center of the Braunbeck coil assembly and whose main purpose was to trigger the 
averager. 

The following conditions were established: 

a. Torquemeter set up (Figure 5) 
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Torquemeter Filter Averager Recorder 

Figure 9. Block Diagram of the Electrical Setup used looking for the 
optimum Operational Frequency. 


b. Flux gate magnetometer near the center of the Braunbeck coil assembly. 

c. Torquemeter load: 400 pounds, 

d. Filter: band pass mode. 

The averager is a sophisticated device whose main purpose is to extract the signal 
from the noise and as a consequence to improve the signal/ noise ratio. This can be done 
when we have a random background noise and the signal is repetitive periodic or not. It 
can also extract the signal from the noise when the period of the noise is different from 
the period of the signal. It has a time disadvantage because you have to wait for a cer- 
tain number of sweeps before the signal can be extracted with some confidence. 

As previously mentioned, the main goal of these tests was to choose the optimum 
frequency of operation. If the natural frequency of operation was chosen, the peak noise 
(See Table 1) was a maximum but the magnification factor is the highest. On the other 
hand operating at 1.00 Hz (Table 1) the noise was minimum but the signal/noise is very 
small. 

A rotational field 0,1 gauss = lOK gamma, ai^ular frequency 0.283 Hz (clockwise) 
and in a horizontal plane was set by the Braunbeck coil assembly and the radial coil was 
energized with a current of 0.640 amp. The magnetic dipole moment of the radial coil at 
this current is: 


M r = 2380 0.640 amp. = 1500 pole-cm 

Amp 


and the expected peak torque is 


T = SM r = 0.1 gaussx 1500 = 150 dyne-cm 


This torque was considered too high. 

Operating the radial coil at 0.064 amp the expected torque is of course t = 15 dyne- 
cm so that the output of the averager should be 0.1 of that obtained when operating at 
0.640 amp. Actually after 20 sweeps the output was found to be 0.08 rather than 0.1. 
Next we operated at a frequency of 1.00 Hz, both the filter and the rotational field and 
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with currents of 0.064 and 0.640 amp in the radial coil. After 80 sweeps the output of the 
averager at 0.064 amps was only about 0.05 instead of the theoretical value of 0.1. There- 
fore the best operational frequency was considered to be the natural frequency (0.283 Hz). 


5. OBJECTIVES AND BRIEF DESCRIPTION OF THE MAGNETIC TESTING 
DONE ON THE SAS-A SPACECRAFT 

Objectives - The SAS-A satellite is the result of the combined effort of the following 
organizations: 

a. AS&E - a research organization in space exploration and located in Boston, 
Massachusetts. Responsible for the X-Ray experiment. 

b. APL - another research organization stressing also on space exploration and lo- 
cated on Route 29, half-way between Washington, D.C. and Baltimore, Md. 

c. The Magnetic Test Site - a branch of NASA/GSFC and whose cooperation in the 
project has been described before. 

d. A private organization located in the Washington area took care of the program- 
ming computer assistance for the computer on board of the SAS-A. 

The main goal of the SAS-A satellite is to conduct an X-Ray experiment looking for 
X-Ray sources coming from celestial sources which apparently are located around the 
equatorial plane. This is the reason for placing the SAS-A in a circular orbit around 
the equator. If the results of the experiment are satisfactory, the SAS-A will be 
followed by the SAS-A but looking for gamma sources. Dr. Riccardo Giacommi from 
AS&E was the principal investigator on the X-Ray experiment. Both SAS-A and SAS-B 
are not intended to be X-Ray or Gamma- Ray telescopes. 

In more detail the scientific objectives of the X-Ray experiment are: 

a. To conduct a high sensitivity, high resolution survey of X-Rays sources, that 
will produce an X-Ray catalogue for sources with intensity greater than approxi- 
mately 5 X 10“* SCOX-1 (the strongest source known at present). 

b. To look for time fluctuations in X-Ray intensities over periods of minutes to 
months. The satellite will be sending information for a period of more than six 
months. 

c. Determination of the spectral distribution of the sources detected in the energy 
range 1-20 KeV. Weight of the experiment is 150 pounds. That of the spacecraft 
is about 180 pounds. So the total weight is approximately 330 pounds. See Figure 
10 for other details of the SAS-A. 

The SAS-A will also keep on board the following facilities: 

a. Star Sensor - consists of a lens that focuses light from the stars on a sensitive 
photomultiplier. Signal output of the photomultiplier will be amplified, trans- 
formed into a pulse and processed into telemetry. 

b. Solar Sensor - The sun sensor will provide information about sunlight when the 
star- sensor is not in operation because of proximity of the satellite to the sun- 
SAS-A line. The sun sensor will also provide attitude information along with 
three magnetometers on board if the star sensor fails to operate. 
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Figure 10. Approximate Longitudinal Crossection of the SAS-A 


c. Star Shutter - The main purpose is to protect the photomultiplier's cathode (at- 
tached to the star sensor) from light coming from the sun when the star sensor 
is pointing towards the sun. 

d. Power System - The power system is mainly provided from solar cells connected 
to the solar paddles. There are 4 solar paddles aligned as shown in Figure 10. 
Only 2 appear in Figure 10. The 2 other solar paddles are also aligned and 
perpendicular to the ones shown in Figure 10. The paddles have the same length 
(approximately 6 feet long) and with their axis in a plane perpendicular to the Z 
axis of the satellite. The cross-section of the SAS-A is circular with its center 
located at the Z axis. 

Comii^ back to the power system it consists, as we said before, mainly of the solar 
array, but there are also a nickel cadmium battery (6 amp-hr) and a shunt regulator. 
Some loads receive power from the main bus and the shunt regulator. Unregulated loads 
receive power through a DC-DC converter. The logic microcircuits of the command and 
telemetry systems are powered through a 5 volt regulator. The command system is also 
provided by another regulator to increase reliability. 

When the power system is operating normally the solar array will be sufficient to 
power the load and charge the battery when it is in the sun light phase of the orbit. 
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The assembly of the satellite was done at APL research center. F. F. Mobley and 
B. E. Tossman from the APL staff were in charge of the stage of the process. They 
were also in charge of the test plan for the SAS-A Pre-Vibration Test. Some innova- 
tions were also made by the personnel of APL. For example, the nutation-damper sys- 
tem, whose main purpose is to contribute to the stabilization of the satellite while flying, 
was designed by B. E. Tossman. J. A. Ford was responsible for the operation of the 
magnetic-trim system related to the attitude control of the SAS-A. To smn up, APL 
took care of and designed the proper hardware for attitude and stabilization controls of 
the satellite. 

In addition to the magnetic testii^ conducted at the Magnetic Test Site (GSFC), the 
satellite will be submitted to Vibration, Thermal and Mechanical Tests. 

The magnetic testing was conducted by personnel of the Magnetic Test Site under 
supervision of J. C. Boyle (GSFC) and personnel of APL xmder direction of F. F. Mobley. 

The testing was done in the 40 Foot Braunbeck Coil Facility. 

A brief description of the Magnetic Testing follows. Data is not included for the 
following reasons: 

a. Data obtained must be reduced, processed and interpreted by personnel of GSFC 
(Magnetic Test Site) and also by APL. 

b. Complete reports on the conclusions of the Magnetic Testing will be prepared by 
personnel of the Magnetic Test Site and APL, but step a, above, must be con- 
ducted first. 

c. The total number of magnetic tests was over 200. Evidently a description of the 
whole process will make this report excessively long. 

For those interested in a complete description of the whole procedure, they can be 
referred to APL Report S2P-2-367, F. F. Mobley and B. E. Tossman. 


MAIN STAGES OF THE PRE-VIBRATION MAGNETIC TESTING 
The whole procedure was divided into four phases: 

Phase 1 - Test-Dipole Moment Measurement by Near- Field Analysis. The goal of 
this phase was to extract the dipole from the other multipoles associated with the 
satellite. Reference is made to report X-325-69-350, W. L. Eichhorn (Magnetic 
Test Site). 

Phase 2 - Test- Dipole Moment on the Torque Table. The main goal of this stage 
was to obtain the dipole-moment associated with the satellite. 

Phase 3 - Test-Magnetometer Alignment, Calibration and Investigation of Biases. 
This phase, that was by far the longest, (2.5 days), was conducted primarily by 
the APL staff with the assistance of the Magnetic Test Site personnel. The main 
goals of this phase are explicitly stated in the title. 

Phase 4 - Spin-Despin System. This took just 0.5 days and the main goals 
were:* 


*F. F. Mobley, B. E. Tossman, APL Report S2P-2-367 
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a. Confirm proper torque levels. 

b. Check insensitivity to bias field. 


GENERAL ADDITIONAL COMMENTS 

a. Launching date: About November 1970. 

b. Launching Site: Launching Platform, San Marco Island, near Kenya, Africa. 

c. Orbit: Circular and Equatorial (300 miles radius). 

d. Expected Life Gathering Data: From 6 months to one year. 

e. Prime Data Gathering and Command Issuii^ Station: STADAN "Space Tracking 
and Data Acquisition Network," Quito, Ecuador, Central America. 

f. Inclination of Orbit: 2.9° 

g. Spin Rate: 1/12 RPM 

h. Period: 96 min. 

i- Time Spend on Preparatory Calibration for the SAS-A Pre-Vibration Test: 6 days. 
]• Time Spent on Pre-Vibration Magnetic Testing of SAS-A Spacecraft: 7 days. 
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STABIUTY OF SAS-A DUAL SPIN SPACECRAFT WITH 
ENERGY DISSIPATION ON THE MOMENTUM WHEEL 


Peter M. Bainum 

Department of Mechanical Engineering 
Howard University 
Washington, D. C. 


ABSTRACT* 

The attitude stability of the SAS-A satellite with damping in the 
momentum wheel as well as the "despun" portion is analyzed. Wheel 
energy dissipation is modeled by assuming the wheel can flex with 
two degrees of freedom relative to the hub. The nonlinear attitude 
equations are derived for small wheel flexural motion and are a ninth 
order nonautonomous set. If the main body damper mass and wheel 
transverse moment of inertia are assumed small when compared 
with main satellite masses and inertias, an averaging process can 
be used to determine the zeroth and first order secular perturba- 
tions on the behavior of the system nutation angle. From this a 
general analytic stability criterion is established. A numerical 
evaluation of this criterion using SAS-A parameters and measured 
wheel damping data indicates that stability about a zero degree nu- 
tation angle is insured by a factor of 128 under normal operating 
conditions. 
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LI, K, RB, SR, BA AND RARE-EARTH CONCENTRATIONS, 
AND RB-SR AGE OF LUNAR ROCK 12013 


Michale L. Bottino 
Department of Geology 
Marshall University 
Huntii^on, West Virginia 


The bulk of this summer's activity evolved working on lunar samples. The investi- 
gation included trace element geochemistry and Rb-Sr geochronology on Apollo 12 soils, 
rocks, and mineral separates. The trace elements studied were Li, K, Rd, Sr, Ba and 
the rare-earth elements. In addition, the Sr isotropic composition was measured. My 
primary area was the Rb-Sr geochronology. 

One paper has been accepted and will be published by the end of August in the journal 
"Earth and Planetary Science Letters." A second, more complete paper will be finished 
by early fall. This work was done jointly with Drs. Charles C. Schnetzler, Paul D. 
Fullagar and John A. Philpotts. The introduction and abstract of the paper are repro- 
duced below. 


INTRODUCTION 

The concentrations of Li, K, Rb, Sr, Ba and nine rare-earth elements (REE) and 
the isotropic composition of Sr have been determined, by mass spectrometric stable 
isotope dilution, in a number of fractions of lunar rock 12013, 10. The primary pur- 
pose of this paper is to present the data obtained to date, and secondarily discuss some 
preliminary conclusions and speculations concerning this rock. 

Examination of this rock in the Lunar Receiving Laboratory showed that it was quite 
high in a number of trace elements, and the mineralogy and general chemical composition 
led the Preliminary Examination Team to the conclusion that it resembled a late-stage 
basaltic differentiate (1). To a first approximation of this rock is a complicated mixture 
of light and dark- colored material. Our first analyses of two prepared, homogeneous 
poweders representing these types indicated that some trace elements had markedly 
different concentrations in the light and dark-colored portions. We therefore decided 
that analyses of the complete small chips we received would not be as instructive as a 
study of separated dark and light-colored portions and some mineral concentrates. 


ABSTRACT 

The light and dark-colored portions of 12013 are characterized by two distinct trace 
element patterns. The light-colored portions are lower in rate earths, slightly lower in 
Sr, the same in Li and higher in K, Rb and Ba when compared to the dark- colored por- 
tions. The dark material appears to represent a late liquid in the normal igeneous dif- 
ferentiation scheme suggested by the other lunar rocks; the light-colored material may 
have had a more complex history. Comparisons of trace element abundances of Apollo 
12 igneous rocks and 12 soil suggest that the dark material of 12013 may constitute an 
important component of the soil. The best-fit isochon to nine data points yields a Rb- 
Sr age of 4.1 + 0.2 b.y. with an initial rate of 0.704 + 0.002. 
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INFERENCES REGARDING THE CHEMICAL EVOLUTION 
OF IGENEOUS LUNAR ROCKS 


John H. Carman 
Geology Department 
The University of Iowa 
Iowa City, Iowa 


A new graphic technique was developed to examine the chemical evolution of 
igneous lunar rocks. The general components for this technique are the metal (M)- 
oxygen (0) complexes MjO, MO, MjOj and MO^. These components represent a 
quaternary system in which the bulk of all rocks and minerals may be plotted. For 
lunar rocks two important exceptions are iron (Fe°) and troilite (FeS). These are 
both treated as the component MO which generalizes it to M, MO, MS. Within this gen- 
eral quaternary the mutually exclusive and exhaustive molecular components 
Kaliophyilite (KAlSiO^), Nepheline (NaAlSiO^) Anorthite (CaAljSi^ 0^ ), the M, MO, MS 
component and the MO^ (SiO^, TiOj, ZrO^ , etc.) component have been employed. Most 
of the compositions examined contained CaO in excess of that required to constitute 
anorthite. This excess, signified as "CaO" was combined with MgO, FeO, NiO, etc. in 
the MO component. Only in the case of lunar glasses (interstitual and spheres) and, 
so called, "anorthosites" was these excess M^O mainly AljOj , over the amount 
necessary to constitute kaliophyilite, nepeline and anorthite. These peraluminous 
compositions are of special importance in regard to the concentration of trace ele- 
ments as indicated below. 

Projections onto the following component planes, Anorthite + M, MO, MS + MO^ , 
Kaliophyilite + Nepheline + M, MO, MS + MO^, TiOj + M, MO, MS + Si02, Nepheline + 

MO^ + Kaliophyilite and FeO + MgO + "CaO" suggest the conclusions that: 

1. The compositional variations of Apollo 11 and 12 igneous rocks (lavas) may be 
explained on the basis of subtractive - and reactive crystallization in which 
minor accumulation of. early crystals or separation of early formed crystals 
are the two most important (and opposing) effects. Accumulation or separation 
of spinel (essentially ulvospinel) and olivine appears to have been important 

in Apollo 11 samples whereas assumulation or separation of troilite and/or 
iron and olivine was most important in Apollo 12 samples. 

2. There is no evidence that there has been any significant early accumulation or 
separation of either plagioclase or ilmenite as has previously been inferred to 
explain the Europium content of these rocks and mascons, respectively. 

3. Maria compositions inferred from these two landing sights are very similar in 
composition and resemble tschermakitic and titaniferous augite cores of Apollo 
11 samples. Available data from experimental studies at high pressures and 
high temperatures suggest that such a pyroxene may be expected to occur in the 
moon's mantle. Thus it seems that such a pyroxenite mantle may be the best 
source material for the maria lavas, assuming complete melting of such material 
during the maria forming stage of lunar evolution. Others have suggested a 
pyroxenite mantle on the basis of density considerations. 

4. The accumulation of troilite and/or iron in Apollo 12 samples suggests that these 
dense minerals may account for the mascons associated with maria centers. It 
cannot be ascertained whether or not the iron and sulfur involved in these miner- 
als is primary or metoritic in origin, however it seems apparent that both are 
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homogeneously distributed within the rocks in which they occur. Significant 
amounts of iron and troilite were also encountered in the Apollo 11 igneous 
rocks, 

5, The K20/Na2 0 ratio has been shown to correlate positively with Th, Rb, La, Ce, 
Ba, Y and Zr, Potassium and these trace elements are known to be concentrated 
in residual liquids (now interstitual glasses) which commonly line vesicules of 
these lunar lavas. Analysis of these glasses indicate a surprisingly low sodium 
content relative to most terrestrial analogs of this type of fractional crystalliza- 
tion in which the glasses approach rhyolite in composition. It is suggested that 
the low sodium contented of these glasses resulted from a fractionation of sodium 
into the gas phase during vesiculation and subsequent separation. This loss of 
sodium enriches these glasses in potassium and the other elements mentioned 
and if loss is extensive it results in excess AljOj in the glass's normative com- 
position, In light of the minerals encountered on the moon there is not other 
known mechanism for obtaining such peraluminous compositions except via the 
fractionation of alkalis into such a exsolvii^ gas phase or the remelting of lunar 
rocks and minerals under the near perfect lunar vacuum. Correlation of 
K 20 /Na 20 with porosity, degree of vesiculation and bulk density, where data are 
available, suggests that it increases in rocks according to their proximity to the 
lunar surface during solidification. These observations and deductions indicate 
that those compositions with low potassium and low trace element abundance 
more nearly represent primary lunar lavas. 
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LOSS TANGENT IN CdTe AT 3,7 GHZ 


Kenneth H. Carpenter 
Physics Department 
East Tennessee State University 
Johnson City, Tennessee 

Frederic M. Davidson 
Electrical Engineering Department 
University of Houston 
Houston, Texas 


ABSTRACT 

An upper limits was found on the loss tangent (imaginary di- 
vided by real part of dielectric coefficient) of CdTe by placing a 
CdTe crystal in a microwave cavity resonate at 3.7 ghz. The value 
of loss tangent was less than 0.2 for a newly prepared crystal, but 
increased to 0.6 after the crystal was subjected to high microwave 
power. The same results were obtained for a second crystal 
specimen. 

A CdTe crystal of dimensions 3x1x1 mm was placed inside a microwave resonant 
cavity in such a manner as to have the microwave electric field directed along the 0,0,1 
direction in the crystal— one of the long faces of the crystal. The long dimension of the 
crystal was the 1,1,0 direction. The entire crystal was situated in the cavity so as to 
experience essentially the peak microwave electric field. The crystal was supported 
in the cavity by being placed inside a dielectric tube of glass or teflon inserted across 
the cavity. The cavity reflection coefficient was measured for frequencies near 3.7 
ghz, the resonant frequency, with and without the crystal inside the cavity (but with the 
support tube always inside the cavity). From these measurements an estimate of the 
dielectric loss tai^ent of CdTe can be made at the 3,7 ghz frequency region. 

The loss tar^ent of a lossy dielectric is defined as follows (see: S. Ramo & J. R. 
Whinnery, Fields and Waves in Modern Radio , 2nd ed., John Wiley & Sons, 1953, p. 

309). For loss expressed as imaginary part of dielectric constant, e = Sg (e' + j e"), the 
loss tangent is e" /e' . Now e" = o-/a> e where o- is volume conductivity, A uniform 
conductor in a uniform electric field of effective (rms) magnitude E has a power dissi- 
pation per unit volume of E ^ ct. If the conducting volume is v , the total power dissipated 
is E^cr V. The cavity Q is given by 


peak energy stored 
energy lost per cycle 


For a rectangular cavity in a 1,1,0 mode we have 


peak energy stored 


- 8 ' 
g 0 max 


v 


where V is cavity volume and is peak value (in time and space) of electric field in 
the cavity. Also S = e' e'E provided we neglect the distortions in the field 
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caused by the crystal and its supporting tube. Thus 
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where Q is cavity quality factor at ai^ular frequency ca with the crystal inside, and Qo 
IS the value with the crystal removed. The value of Q was measured by observing v^i- 

ation of reflation coefficient with frequency and hence finding the half power frequencies 
andthereby QfromQ = f/Af. We thus find luquencies 


= i [(Af) 

e 4 V f 


(A f)J 


A - essentially identical CdTe crystals were supplied by James Keifer of Hughes 

Aircraft Company. The toss tangent was measured for one of the crystals in the condi- 
tion in which it was received. Both crystals were used in a laser modulation experiment 
during which they were subjected to microwave powers on the order of one to less than 
en watts Durii^ this time the toss tangent increased significantly. Measurement of 
the toss tangent for each of the crystals following the appUcation of high power showed 
essentially the same toss tangent, which was at least three times the initial value. 

formula for toss tangent to the crystals as received is difficult 
since loss in the tube used for support was much greater than in the crystal. Hence 
change in A f was less than the accuracy with which A f was determined. Assuming a ten 
percent accuracy in Af, we can bound the loss tanget in the new crystal as 


< 0.1 i XU 


4 v/Q 


The numerical values are 


V = 1 X 1 X 3 X ( 10 ' 3)3 m 3 


V = 2.5x5x5x (10-2)3 m3 
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Hence 


i X = 5.2 X 1Q3 
4 V 


0= 2x 1Q3 


— < 0.1 X 
6 


5x 103 
2 x 103 


0 . 2 . 


After the crystal had been irreversibly damaged when between one and ten watts of 
power was incident on it, the losses were increased and measurements gave Ff^ = 1.7 
mhz andAf^j^i = 2.1 mhz at f = 3.7 ghz. Thus 


e" 

7 ^ 


^ [2.1 - 1.7] X 10« 

3. 7 X 109 


0.6 


Thus the loss tangent had been increased by at least a factor of three, and probably 
much more, by whatever mechanism resulted when high fields and, or temperatures were 
reached by the crystals. 
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ABSTRACT 

Mode coupling in an optical resonant cavity is examined from the 
viewpoint of modulation of a single mode laser. It is shown that an 
intercavity modulator can achieve enchancement of the modulation in- 
dex over external modulation by a factor on the order of the passive Q 
of the cavity. Such internal modulation may be achieved over bands- 
widths exceeding that of a single resonator mode by coding the modu- 
lating singal into frequency bands about successively higher resonant 
modes. 


Amplitude modulation of laser beams through the use of external modulators has 
been reported at various frequencies up to several gigahertz Internal modu- 
lation of single mode lasers, either AM or FM for information coding purposes, has 
generally been limited to frequencies less than the optical cavity line width of a few 
hundred megacycles. The use of internal modulation techniques is preferred because 
a lower modulator drive power is required for a given power in the sidebands. Hence 
it is desirable to extend internal modulation techniques to larger bandwidths than that 
of a single optical cavity mode. In the following we will consider internal modulation 
in an optical resonator from a mode couplir^ viewpoint, and will derive an expression 
for the modulation index which explicitly shows modulation is enhanced over external 
modulation by the Q of the resonator. 

GENERAL EQUATIONS FOR INTERNAL LASER MODULATION 

Internal modulation of a laser can be understood from the viewpoint of mode coupling 
in a resonator. This problem has been extensively studied by Yariv and others. ^ In 
the following we use the method and notation of Yariv®, corrected and adapted as neces- 
sary to fit our application. 

First, we define for the optical resonator (cavity) the orthonomal sets of vector 
functions (7) and (T) satisfying 


k E = Vx H . k H = Vx 

a Q a a a a 


( 1 ) 


V-E„ = V-H = 0 

0 a 


(2) 
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and the boundary condition n x =0. Orthonormality implies 


f E,-E,dv= [ H^-H, 

♦'V *'v 


d V = S 


(3) 


Either set or set is complete for expansion of a vector field in the cavity that 
satisfies the same boundary condition as or . Now for an optical cavity with walls 
of perfect conductors the electric field vector satisfies the n x e = 0 boundary condition. 
We will assume this is the case to a good approximation and account for actual losses in 
the walls by a fictitious contribution to volume conductivity. Hence, we expand the 
electromagnetic fields inside the optical resonant cavity as 


E(r, t) 




Pa ft) K 


(4) 


H (T, 


t) . 


"a 


(5) 


where e and are the effective homogeneous constant total permeability and permittivity 
with the modulator present in the cavity but with no drive to the modulator.* By definition, 
“a example, in a Fabry-Perot interferometer of length L, = 7 rn^/L. 

Next we use Maxwell's equations 


V X E = - , 


3H 
3 t 


( 6 ) 


Vx Hr J + 4 [e(7. t)f (7, t)] 

0 t 


(7) 


with expansions (4) and (5) to show 


and 


Pa(t)=^ [q, (t)] 




"S ‘lb + ■ 


t) 


Pk + 


3 t 


£ (r, t) p. 


\(T)=0 


( 8 ) 


(9) 


When the modulator does not fill the entire cavity the static e is not homogeneous. However, the net 
effect is that of a homogeneous £ of appropriate value since static inhomogeneities do not couple 
modes. 
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where we have expanded current density as 


Now 


J = 



Pa *^a 

JT 


(10) 


CTfa ( r, t) 
€ 


"b <^ib 

" Qb" ^ 


( 11 ) 


with Qb being the total effective Q of the resonator in mode b in the absence of modulator 
drive (including losses due to finite conductivity walls, output coupling, and constant vol- 
ume conductivity) and a^ being the time varying conductivity induced by modulator drive. 
Also we define (7 , t) by 


e (r, t) = e + ( r. t). 


( 12 ) 


The b subscript is placed on and in (11) and (12) to indicate that if the modulator 
medium is dependent on optical frequency, then in each mode of the expansion in (9), Ej or 
cTj must be evaluated at the actual frequency present in the mode bandwidth. Thus on 
takii^ a scalar product of Eq. (9) with and integrating over the cavity volume we find 


where 


and 


CO 

Pa + -^ Po + "a ‘la 



s 


ab 



t) E • E. d V 

' a D 


d , 1 

d t e 



t) E • E. d V 

' a b 


(13) 


(14) 


(15)* 


Equations (8), (13), (14) and (15) form the general basis for mode coupling in a 
resonator due to internal modulation. We next specialize to the case of sinusoidal modu- 
lation of a single mode laser and proceed to solve the resulting equations. In order to 
obtain the solution the following assumptions are made: 


*The effect of loss in o dielectric at optical frequencies is often stated by giving an imaginary part to 
s (T,t). This follow's from an assumption of a single frequency complex form forE (r,t) = C ei“‘ and 
from the time derivative operator in Maxv^ell's equations yielding a factor jcj. Thus when a complex 
[e (T, t) is given one must use Re [Ej( r , t) for €)(■■/ t) in Eq. (14) and use m L£j(r , t) for 
o'i( T',!;) in (15). 


33 



(16) 


1. S^^ = Slsinvt = -^j (eivt_e-jvt) 

2, (t) = cos (v t + Q') = (17) 


3. The single mode laser oscillation is on mode a = L and has constant /. 

amplitude, (t) = A cos t 

4. The upper sideband frequency (oj^ +v) lies within the pass band of 

mode a = u, i.e., + v falls in (1 ± 1/q^ ), and all other possible 

frequency components, e.g,, - v , o) ± 2 v , etc, lie outside 

(1 ± 1/Q„). The corresponding assumption is made for the lower 
sideband in mode a = -t. 

Using these assumptions we proceed to specialize the differential equation for the 
upper sideband mode. First we eliminate from (13) with (8) and then write (8) for 
a = u to obtain 


+ ^ q + co; 




+ s 


3, + R 

uu 


% + Ki 


(19) 


With the assumed values of q, , R . , S ^ we have 

do do 




CO 

+ 5 -” ‘5u + "u % = - ^ "l RSl cos (V t + 0) sin t 
- A CO? S° sinvt COSO), t 

L uL 2 

+ terms involving sin v t or cos v t times sin (co^t v) t or cos co^ 


( 20 ) 
t vl t 


By assumption 4 above the only driving terms on the right of (20) that need be retained 
are those of frequency (<^l +'')• Hence we may approximate (20) by 


CO 

.. u • 9 

Q -I- ' — q + q 

^ 0 u 

( 21 ) 

= - i A sin [(ojj^ + v)] + R®|^ sin i(co^ + v) t + Q] } 


Equation (21) is the standard driven harmonic oscillator with the solution 


A G sin + v) t + <fi 


[(co^ + v)2 - oj2]2 + 


(co + v)2 oj2 


q„ (tl = 


1/2 


( 22 ) 



where 


\ 


\ 

■ \ 
I 


and 


G . 1 + 2 ..j, SO^, R«, cos 4>y 


1/2 


(23) 


Kp = tan” 


/ sin \ _ / ("l + ''^ "u 

tan”l / + tan > 


R^l <^os4>+co^ 


(K + v]2 _ 0)2-) 


\ (24) 


Details of this solution are given in Appendix A. 

The ratio of electric or magnetic field strength magnitude at the upper sideband to, 
that of the carrier (laser oscillation) may now be calculated. From (5) and (22) this \ 
ratio is 


(25)\ 


Using G from (23) we find 


M , = G J [(O), + v) 


(o)^ + v ")2 0)2 


- 1/2 



The value of M in (26) gives the ratio of field strengths inside the laser cavity. To 
find the ratio outside we must multiply by the ratio of output coupling for the two modes. 
In (26) Qu is the total effective Q for mode u including the effect of losses in the modu- 
lator. ’Thus the ratio of field strengths outside the cavity will not be QyQ^ because 
the Q and Q, values include losses other than coupling. Instead we cau write for the 
ratio of field strengths outside the laser m^^^ = (tu^L^ where t, is the transmission 
coefficient of the output mirror at the appropriate frequency. 

Essentially identical results will be obtained for the lower sideband by substituting 
the appropriate values into Eq. (26). 


DISCUSSION OF INTERNAL MODULATION 

The expression for resultant modulation index for internal modulation given in Eq. 
(26) shows the advantages to be achieved by this method of modulation. First, we should 
note that at the center of the cavity mode, >_the ratio of power in the upper 

sideband to power in the carrier in the external beam is 





R“, 


+ 2 — cos 4> 


(27) 
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The significance of this ratio is in the presence of the factor Q 2 which can be quite large 
for an optical resonator. However, for loss modulation, ^ (7, t) % 0, the effective Q could 
be small if there are a D.C. (unmodulated) component of ct( 7 , t) at the uppersideband 
frequency. Hence, if the enhancement of modulation due to presence of the resonator is 
to occur with loss modulation the loss mechanism must be frequency dependent and losses 
must be confirmed to a band containing but not containing + v. In fact, when the 
modulation index is calculated assuming a frequency independent loss mechanism the 
factor of Q cancels in the result. From this argument we see that the factor in 
(27) enhances the sideband power significantly (i.e., by a factor of 10^ to^’lD^o ) over that 
obtained by external modulation with the same modulator and modulator drive provided 
the modulator losses are low at the sideband frequencies. The method should therefore 
be useful for both Stark effect loss modulation, s. and for electro-optic modulation of 
e (r ,t) when the electro-optic modulator is oriented so as not to produce losses due to 
polarization effects. 

Although internal modulation enhances the modulation index by the resonator Q when 
the sideband frequency falls on a resonator mode center frequency, for other sideband 
frequencies the modulation is reduced by the factors in the denominator of Eq. (26) in 
the manner of a second order resonant circuit with quality factor and center frequency 
• Thus in a given mode of the resonator the bandwidth of modulation frequencies that 
can be used is and the enhancement of the gain-bandwidth product for the modula- 

tion is Qu(“u/^u) ~ "u is independent of Q^. Hence in this sense no enhancement is 
achieved by internal modulation over external. But, with high Q, the bandwidth limitation 
on internal modulation may be circumvented by coding of information to place the modula- 
tion frequencies in successive bands of width co/Q at successive resonator modes, thereby 
achieving both unlimited bandwidth and the enhancement of modulation by Q. Thus in- 
ternal modulation appears attractive for gigahertz frequencies where modulators are in- 
efficient without such an enhancement. 

One may foresee problems with maintaining stable laser oscillations with an inter- 
cavity modulator, particularly a loss modulator. However, the derivations above will 
apply equally as well to the case where the modulator is placed in a separate optical 
resonator from the laser resonator thereby isolating the modulator from the active 
medium. (This would also be necessary if the active medium has gain at sideband 
frequencies if distortion is to be avoided.) 

In summary, one may achieve an enhancement of laser amplitude modulation by the 
Q factor of an optical resonator by placii^ the modulator inside the resonator (provided 
the resonator Q at the sideband frequencies is not significantly degraded by the modulator). 
Such internal modulation may be achieved over bandwidths exceeding that of a single 
resonator mode by coding the modulating signal into frequency bands about successively 
higher resonator modes. Thus it should be possible to achieve gigahertz and higher 
modulation of lasers with internal modulators. 
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APPENDIX A 



(oi^ - (A sin 6j' t + B cos oj' t) + C- B sin oi' t + A cos oi' t) = a sin co' t + /3 cos a>' t 

A 

which on equating coefficients of cos a>' t and sin o>' t yields 



or on inverting the matrix 
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Next we substitute the appropriate values for w , co' ,a , p for Eq. (21) to find G, i/<of (23) 
(24): 

a = - 1 A cos 0] yS = - i A sin 0 

a2 + ^2 ^ ^ A2 0)2 {0)2 ^ ^ 2 R^l (A5) 

B _ _ - ("l + ''■> "u K ^uL + I^L + ^v, K - C"l + ’^uL 

^ K - ("2 + K ^UL + ^UL COS + (o)j^ + V-) O)^ ROj_ sin 4y 

Since tan ^ = B/A we may write 

R®l sin </> (o)j_ + v) o)^ 

"2 s!l + Kl ^ Qu [("l + v)^ - ■ 

tani/-= (A6) 

SSl + R®l cos ./.j (p^ [( 0)2 + v)2 - o)2]j 

and use the trigonometric identity 

tan (A + B) = tan A + tan B 

1 - tan A + tan B 

to obtain (24) from (A6), Eqs. (22) and (23) follow at once from (A4) and (A5). 
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EVALUATION OF COMPUTER STRUCTURES FOR THE STADAN 
AUTOMATED DATA HANDLING AND CONTROL SYSTEM 


Stephen W. Ching 

Department of Electrical Engineerii^ 
Villanova University 
Villanova, Pennsylvania 


This summer I am working in the evaluation of computer structures for the STADAN 
automated data handling and control system. My GSFC associate is Mr. John C. Rockers. 

The STADAN automated data handling and control (STADAC) system is an automated 
system consistii^ of special pxirpose hardware, computers and system softwares planned 
for more effectively perform data handling, station operation and equipment control at 
the remote station of the world-wide space tracking and data acquisition network (STA- 
DAN). The need for evaluation of computer structures arises initially when the need for 
a computer facility is determined. The need for evaluation is never satisfied completely 
thereafter. The plans for implementing a computer facility at the remote station of a 
space tracking and data acquisition network involve the following basic question: "What 
kind of computer configuration is required to perform the anticipated data handling and 
system control task within a required response time?" It is clear that many different 
computer configfurations could satisfy the STADAC requirements. The objective then, 
is to determine which computer configuration is optimal. The optimal computer structure 
of STADAC must be considered relative to the six STADAC functional subsystems such 
as data handling subsystem (DHS), schedule display subsystem (SDS), equipment status 
reporting subsystem (ESRS), equipment set-up and control subsystem (ESCS), link control 
monitor subsystem (LCMS) and link readiness and verification subsystem (LRVS). In 
order to make a meaningful evaluation of various computer configm-ations, standard 
measures of system capabilities must be employed. Some of the measures related to 
the STADAC requirements are turn- around- time, throughout, cost, system reliability 
and the combination of these factors. Throughout of STADAC system is measured in 
time units as the capability of the system to accept data from orbiting spacecrafts and 
transmittii^ the data back to GSFC. 

A detailed study of STADAC functional subsystems revealed the following computer 
structure characteristics; 

1. The major computer design crtieria will be optimal (effective) use of available 
communication and control interfaces. The system will be classified as a com- 
munication and control computer and will be capable of handling mass data base. 

2. The system will be controlled primarily by data rather than by program. 

3. Use of hardware to govern communication and control interface will be 
emphasized; microprogramming facility is highly recommended. 

4. Most processing will be executed in real time. 

5. The system will be readily expandable. Hardware and software will be modular 
in design. Computing power will be modified without redesign of the system. 

6. The operating system will be designed to operate efficiently. The system will 
include an analysis of I/O buffering requirement, page characteristics or mem- 
ory fragmentation, time slicii^ algorithms for multiprogramming, queuing 
disciplines as applied to scheduling and dynamic aUocation of system resources. 
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Three kinds of computer configuration are considered in this evaluation. The 
first configuration is a dual processors multi- link computer system which contains a 
data-handling central processor unit and a station control CPU. A second configuration 
is a link-oriented design using mini- computers. Array computer with micro-facility is 
then considered. Evaluation results of the above computer configurations will appear 
in my final technical report. 
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I. INTRODUCTION 

Modulation of a CW CO, laser at gigacycle modulating frequencies was attempted 
in two ways. The first was through the use of a crystal that exhibits a linear electro- 
optic effective, i.e., on that has an index of refraction that can be varied linearly with 
applied electric field strength. The second method involved the use of a gas cell modu- 
lator containing a gas that exhibits a Stark Effect energy level splitting for energy levels 
that give rise to absorption lines in the 10.6 micrometer band of CO, laser transitions. 
Theoretical considerations for both types of modulation will not be considered here, but 
a general theoretical discussion of laser modulation techniques can be found in the sec- 
tion of this report entitled "Internal Modulation of Single Mode Lasers at Frequencies 
Exceeding Optical Cavity Line-Width." 


II. LINEAR ELECTRO-OPTIC EFFECT MODULATOR 

Modulation of CW visible laser beams at gigacycle frequencies has been demon- 
strated with the use of Li NbO, and also with KDP^. Modulation of a CW CO 2 laser at 
10.6 micrometers has been accomplished with the use of GaAs*. It was the purpose of 
this investigation to attempt modulation of 10.6 micrometer radiation with the use of 
CdTe since it has a larger electro-optic coefficient than GaAs and hence would require 
less modulator power to achieve the same depth of modulation^. 

Two crystals of CdTe, approximately 1 mm x 1 mm x 3 mm, were obtained from 
the Hughes Aircraft Corporation. CdTe belongs to the 43 m crystal symmetry class 
and is a cubic crystal. The linear variations of with applied electric field is 

described by the electrooptic tensor^ where i=l, 2, 3, 4, 5, 6, and j = 1, 2, 3. 


A 



r.. Ej 
‘J i 


Repeated indices are to be summed over are respectively 
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The coefficients of the yz, xz, xy of the index ellipsoid written in a nonprincipal co- 
ordinate system. Since the crystal is cubic, the coefficients of the X*, Y*, and terms 
of the index ellipsoid are equal and may be denoted by 1/v^ ■ The electrooptic tensor for 
43 m symmetry class crystals is of the form 

/ 0 0 0 r^, 0 0 \ 

41 

0 0 0 0 r^, 0 

41 

0 0 0 0 0 r^, 

\ / 

CdTe has a value of '41 = 5.7 x 10"*^ meters/volt. 

The crystals were cut so that their long axis corresponded to the [110] crystal di- 
rection, and two of the long parallel faces corresponded to planes perpendicular to the 
[001] crystal direction. The orientation of the optical field electric vector with respect 
to the modulating r-f electric field vector is determined by finding in which direction the 
maximum change of index of refraction with applied electric field occurs. This is found 
by finding the intersection ellipse between the plane through the origin perpendicular to 
the direction of propagation of the optical field and the index ellipsoid. To obtain maxi- 
mum modulation the optical field should propagate along the long crystal axis. The r-f 
modulating field can conveniently be applied along the z direction [001] crystal direction. 
The changes of the index of refraction with applied electric field are consquently de- 
termined by the intersection of the plane 


X + Y = 0 


and the index ellipsoid 


+ y2 + z2 


^0 


+ 2 r E X y 


1 


The intersection ellipse is given by 


2 



(1 - r,, E. 



or by 


E. vS) + ZVr,* = 1 

Hence modulation results only for the X and Y components of the optical electric 
field, and the change in the index of refraction in the X or Y direction is given by 


A 77 = 



■■41 E, 
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The phase change of the X and Y components of the optical field produced by a crystal of 
length L are given by 


2 TT 

T 


L A 77 = 


2 \ 




E 

Z 


If the modulating field, Ez is of the form Ez = E^^ sin the X and Y components of the 
incident optical field are phase modulated and the resultant optical field components are 
given by 


E =K sin(&; t + ^ ^ 

Opt L>o \ opt 2 A. ^ 


Ui ^^0 "m t 


Elo is the maximum aplitude of the optical field electric vector after passing through the 
crystal of length L. The ratio of the power in the first side band of the optical field to 
the power in the optical carrier is given by 

k is the wavelength of the incident optical field, and vi '"41 = 10 meter s/volt for 

CdTe. 


The apparatus used to modulate a COj laser with the CdTe crystals is shown in 
Figure 1. The laser used was a frequency stabilized Honeywell flowing gas system. The 
CdTe crystal was placed inside a r-f cavity, resonant at 3.77 GHz with a loaded Q of ap- 
proximately 2000. The modulating electric field was applied through the use of a stable 
r-f oscillator and traveling wave tube amplifier capable of supplying 10 watts to the 
resonant cavity. The frequency components of the modulated laser beam were determined 
with the use of a scanning Fabry- Perot interferometer and gold doped germanium photo 
conductive detector operated at 77°K. 



Figure 1. Apparatus used to phase modulate a CW CO 2 
using CdTe electrooptic modulator. 


laser beam 
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The electric field, E , that could be applied to the CdTe crystal is given by 

*o 


->0(1 -Pr) PjN ’ 

_ fo ^0 V _ 


where 

P.^ = power incident on the r-f resonant cavity 
= power reflection coefficient 
Q = Cavity quality factor 
tp =8.85x10*'^ farads/meter 
fO = cavity resonant frequency 
V = cavity volume 

The power reflection coefficient and cavity Q were determined with the tuning curve 
shown in Figure 2, and a spectrum analyzer. These measurements yielded a value of 
given by E,o = 360 x p, volts/cm where Pi„ is the r-f power in watts incident on the 
resonant cavity. The maximum expected value of the ratio of the power in the first side- 
band to power in the carrier obtainable with the experimental apparatus is given by 


[2 . 1 , ( 0 , 016 ) ( 0 . 016 )] 2 ^ 2.6 x 10 “* 



Figure 2. RF resonant cavity tuning curve. The 
cavity contained a CdTe crystal that had not 
been damaged by high RF povi^er levels. Hori- 
zontal scale is approximately 1.5 MHz per 
centimeter. 
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As r-f power was applied to the crystal of CdTe, both crystals were observed to be- 
come more lossy after the absorbed r-f power levels exceeded 1 watt. The loaded cavity 
Q was observed to change slightly (by less than 10%), but more importantly, the power re- 
flection coefficient increased from 0.1 to 0.8, as shown by the tuning curve in Figure 3. 
This substantially reduced the maximum value of that could be applied to the crystal. 



Figure 3. RF resonant- cavity tuning curve after 
the CdTe crystal dielectric properties had been 
altered by RF power levels in excess of 1 watt. 

Horizontal scale is approximotely 1.5 MHz per 
centimeter. 

The crystals were also observed to absorb heat more rapidly after the applied r-f 
power levels had exceeded one watt. The crystal ends were parallel, polished surfaces, 
and hence the crystal acted as an etalon. The thermal expansion of the crystal caused the 
transmitted optical field to undergo variations in intensity that were easily observable. 

The rate at which these variations occurred increased by more than a factor of 10 after 
the change of cavity Q and power reflection coefficient had been observed. The change in 
the crystals' dielectric properties was irreversible. 

Figure 4 indicates the output form the scanning Fabry-Perot interferometer. The 
peak corresponds to the laser carrier. Sidebands, if present, would be indicated by small 
peaks somewhere on the trace. The translateable mirror in the interferometer was at- 
tached to a piezoelectric cylinder whose maximum extension with applied voltage was only 
1 micrometer. The entire translateable mirror assembly could be moved with a differen- 
tial screw, allowing the interferometer to be manually scanned through one complete free 
spectral range (5.3 microns, or 600 MHz). No evidence of sidebands was ever seen. This 
was due to the extremely low values of the ratio of power in the sidebands to power in the 
carrier and the inability of the CdTe crystals to withstand applied electric field strengths 
of sufficient magnitude to produce sufficient, phase modulation in order to create detectable 
power in the sidebands. 

The optical insertion loss was measured for both crystals, and found to be approxi- 
mately 30%. This insertion loss was unchanged after the crystals were subjected to r-f 
power levels in excess of one watt. The high optical insertion loss precluded the possi- 
bility of modulation with the CdTe crystal inside the laser resonant cavity. 


47 





Figure 4. Ou^put of scanning Fabry Perot inter- 
ferometer. The large peak corresponds to the 
laser carrier frequency. Vertical scale is 1 mv 
per centimeter. Horizontal scale isapproximately 
12 MHz per centimeter. 

in. GAS CELL MODULATORS 

Modulation of CW CO^ lasers oscillating at wavelengths of 9.6 and 10.6 micrometers 
at megacycle modulation frequencies has been observed with the use of low pressure gas 
cells® , Modulation is achieved by variations in time of the absorption coefficient of the gas 
caused by an applied, time varying electric field’ * . If the gas modulator is not located in- 
side an optical resonant cavity and the laser beam simply passes through it, the ratio of 
the power in the first sideband to the power in the fundamental is given by® 

L2 T, 5/,To (o-'M]"- 


where 


a:' = maximum Stark Effect energy level splitting 
j/ = modulation frequency 
T = mean collision time of gas molecules 

s; = frequency difference between optical radiation frequency and dc Stark Effect 
absorption frequency. 

The Stark Effect energy level splitting is given by /j.EM/t) J (J + i) where ,, is the dipole 
moment of the gas molecule for a particular electric dipole transition, E the applied elec- 
tric field, h is Planck's Constant, J the total angular momentum of the gas molecule and 
M the Z component of J. In the event that ^ = 0, no odd harmonic sidebands appear. The 
ratio of the power in the first even harmonic sideband (at frequency 2i/) to the power in 
the fundamental is given by 


[2 T, (a' (cu' 'vS\\ 
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If the gas modulator is located inside an optical resonant cavity, the ratio of power 
in the sidebands to power in the carrier may be enhanced by a factor 




2 


(see the section of this report dealing with internal modulation of single mode lasers). 

Qsb is the Q of the optical resonator at the sideband frequency and tgB , tj, the trans- 
missivity of the resonator mirrors at the sideband and laser fundamental wavelengths, 
respectively. 

The gas cell used in the experiments consisted of a RF cavity resonant at approxi- 
mately 3.8 GHz that had two salt windows set at the Brewster angle for 10.6 micrometer 
radiation attached to it. The cavity could be evacuated and filled with an absorbing gas 
at low pressure. The gas cell was of sufficient dimensions that it could be fitted inside 
a laser resonant cavity. 

The apparatus used to modulate a CO^ laser with the gas cell inside the laser opti- 
cal resonant cavity is shown in Figure 5. The laser oscillation frequency could be adjusted 
through the use of a mirror attached to a piezoelectric cylinder (PZT). Total available 
mirror motion was 12 micrometers. Spectral analysis of the optical field with a scanning 
Fabry- Perot interferometer was again used to detect the presence of modulation of the 
laser beam. The gas used to modulate the laser was 1,1, difluor ethane, C 2 H 4 F 2 , which 
has two very strong absorption lines in the 10.6, micrometer CO, laser transition band, 
one at 944 cm"‘ (P(20)) and one at 939 crn'i (P(26)). 



CO2 LASER 


Figures, Apparatus used to modulate CO2 laser with a gas modulator cel I 
located inside the laser optical resonator. 


An investigation was carried out to determine the effects on the COj laser signature 
caused by the presence of the gas cell inside the laser resonator. Figure 6 gives the laser 
output power as a function of the position of the moveable mirror with no C 2 H 4 F 2 in the 
gas cell. The portion of the signature contained in the fifth box from the left and at the 
extreme right of the trace corresponds to the P(20) laser transition. 
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Figure 6. Laser signature for the case of an eva- 
cuated gas modulator cell. 



Figure 7. Laser signature when the gas modulator 
cell contains Cj F^ at a pressure of 0.2 mm Hg. 

Figvire 7 shows the effect on the laser signature of the presence of Cj at a 

pressure of 0.2 mm Hg in the gas absorption cell. The presence of the gas causes the 
laser to oscillate on frequencies that are not absorbed by the gas, or to oscillate on a 
frequency that is absorbed by the gas but with a reduced amplitude of oscillation. The 
laser will, of course, oscillate at the frequency that has the highest gain in the laser 
resonant cavity. This is determined by the exact length of the cavity, the doppler pro- 
files of the laser transitions, and the absorption losses of the gas modulator cell. The 
P(20) laser transition corresponds to the portion of the signature contained in the fourth, 
fifth and tenth boxes from the left. As the H^Fj gas pressure in the modulator cell 
was increased, the laser output power at the P(20) oscillation frequency was reduced, as 
shown in Figures 8 and 9 (CjH^Fj pressures of 0.3 mm Hg and 0.7 mm Hg respectively). 

A C 2 H 4 F 2 pressure of 1.2 mm Hg was sufficient to cause the laser to stop oscillating on 
the P(20) transition altogether, as indicated by the laser signature shown in Figure 10. 

The laser output at frequencies absorbed by the Cj H 4 F 2 was observed to be ampli- 
tude modulated at a low rate (10 - 100 KHz) as indicated by the signature shown in Figure 
11. No externally applied electric field was applied to the gas in the modulator cell to cause 
this amplitude modulation. The pertinent features of this effect are summarized below. 


50 





Figure 8. Laser signature when the gas modulator 
cell contains Cj Fj at a pressure of 0.4 mm Hg. 



Figure 9. Laser signature when the gas modulator 
cell contains CjH^ Fj at a pressure of 0.7 mm Hg. 



Figure 10. Laser sigrrature when the gas modulator 
cell contains H^Fj at a pressure of 1 .2 mm Hg. 
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Figure 11. Lasersignature showing kilocycle am- 
plitude modulation of laser output frequencies 
that are absorbed by the gas inthe modulator cell. 
The cell contained C 2 F 2 at a pressure of 0.2 
mm Hg. 


1. The amplitude modulation occurs only at frequencies absorbed by the gas in 
the modulator cell. 

2. The frequency and wave shape of the amplitude modulation were extremely 
sensitive to the exact length of the laser resonant cavity. A change of a ap- 
proximately 0.1 micrometers in the length of the laser resonant cavity caused 
the modulation frequency to change by a factor of 10, and the modulation wave 
shape to go from triangular to sinusoidal. 

3. The modulation depends on the absorption line width of the gas. If the absorp- 
tion lines are pressure broadened by mixing the gas with air, the amplitude 
modulation does not occur. Figure 12 gives the laser signature with the modu- 
lator cell containing Fj at a partial pressure of 0.1 mm Hg and air at a 
partial pressure of 10 mm Hg. The signature is nearly identical to the signa- 
ture shown in Figure 8, and differs only in the lack of kilohertz amplitude modu- 
lation of the laser oscillation frequencies absorbed by the H^F^ . 

No satisfactory explanation for this low frequency amplitude modulation was found. It is 
not likely that it is due to the laser switchii^ between different oscillation frequencies, 
since once it switched to an oscillation frequency for which the cavity gain was increased, 
it would remain oscillating at that frequency. The modulation is not due to a laser oscil- 
lation line shape containing two dominant frequency components separated by less than 
the laser cavity linewidth. Such an output line shape would require that the absorbing gas 
have an absorption line width of the same order of magnitude as the observed amplitude 
modulation frequency. The doppler broadened 944 cm"> absorption line linewidth of 
CjH, F 2 at 300°K is 15 MHz. 

Attempts were made to amplitude modulate the laser oscillating at the P(20) oscil- 
lation frequency with the gas cell inside the laser resonant cavity. The cell contained air 
at a partial pressure of 15 mm Hg and CjH^Fj at a partial pressure of 0.2 mm Hg. This 
mixture could withstand an incident r-f power level of about one watt before the gas 
molecules became ionized. This corresponded to an applied electric field strength of 
about 300 volts/centimeter. The resultant Stark energy level splitting, a' was of the 
order 10^ Hz. The modulation frequency, v , was 3.8 x 10® Hz. The expected ratio of 
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Figure 12. Laser signature when gas modulator 
cell contained CjH^Fj at a partial pressure of 
0,1 mm Hg and air at a partial pressure of 10 mm 
Hg. 

power in the first sideband to power in the fundamental was of the order if t x 10'®. 
Even though no dc electric field could be applied SJ could be varied over several mega- 
hertz by tuning the CO 2 laser over the doppler broadened P(20) laser transition. The 
transmissivity of the resonator mirrors was the same for both the sideband and laser 
fundamental frequencies. The optical resonator Q for the sideband, Qgg, may have been 
low since the sideband frequency may not have been close to multiple of the axial mode 
spacing of the laser resonator. 

The scanning Fabry Perot interferometer was used for spectral analysis of the 
modulated laser beam. No evidence of modulation was ever seen. This was probably 
due to the condition a t << 10®. This factor could not be increased since neither 
the modulating frequency u nor the length of the laser resonator cavity could be varied 
so that the sideband frequencies coincided with the multiples of the laser resonator axial 
mode spacing. If the sideband frequencies did coincide with a resonant frequency of the 
optical cavity, Q|g could easily be of the order 10*° to 10 


rv. CONCLUSIONS 

Although neither attempt at modulation of a CO^ laser at gigahertz modulation 
frequencies proved successful, both methods could be made to work. In the case of 
modulation with CdTe crystals, more work is needed to determine techniques for pro- 
ducing CdTe crystals that can withstand high r-f power levels. 

Gas cell modulators have been demonstrated to effectively modulate CO^ lasers 
at frequencies for which /v ^ 1. Modulation at frequencies for which u,‘ / v < < 1 is im- 
practical for modulator cells outside an optical resonant cavity due to the extremely low 
ratio of power in the sidebands to power in the laser carrier. Location of the gas modu- 
lator inside the laser resonant cavity drastically alters the laser signature and could 
prevent frequency stabilization of the laser with the use of the power profile feedback 
type stabilization schemes. However, a gas cell modulator located in a properly designed 
passive optical resonant cavity could prove to be an efficient method of modulation at 
multiples of the optical cavity axial model spacing. 
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SCATTERING AND PROPAGATION 


Garabet J. Gabriel 
Electrical Engineering Department 
University of Notre Dame 
Notre Dame, Indiana 


1. RESEARCH 

The research effort was addressed to the problem of scattering of electromagnetic 
waves from earth and its effects in multipath communication between two satellites such 
as in the proposed Data Relay Satellite System and Navigational Satellite. To date no 
rigorous or semirigorous analyses of the problem, appropriate to conditions of DRS sys- 
tem, exist in the literature. The need for a fresh approach became apparent. Principal 
conclusions of this attempt are summarized next; details will appear in a GSFC document 

(a) The method of current distribution, heretofore considered an approximation 
which is strictly valid only on infinite flat surface illuminated by a plane wave, 
is proved to be one of the weakest foundations of treatments of scatterii^ from 
the real earth. 

(b) When a large, conducting, irregular object is illuminated with a plane wave, 
only the point on the surface which is directly below the receiver contributes 
significantly to the scattered field at the receiver. The asymptotic method 
used here is not conclusive at this stage when the illumination is from a finite 
source. 


2. GSFC COMPUTER APPROACH 

An existing program at GSFC was being adapted for numerical integration of the 
rigorous integrals for it was found, the field oscillates erratically as the integration area 
is increased, no criteria for truncation could be determined. Requested to resolve this 
difficulty, I demonstrated analytically that the scattered field should oscillate with in- 
creasing area, and that integration must extend over entire illuminated region. Accord- 
ingly, it was concluded that it would be impractical and uneconomical from the standpoint 
of time and cost to integrate numerically over the entire illuminated region. 


3. CONSULTATION 

I was invited as consultant-observer to a conference with representatives of GSFC, 
TSC, and FAA to resolve differences in treatments of data on polar cap absorption, scintil- 
lation absorption, and multipath effects as pertinent to the proposed Aircraft Navigational 
Satellite. Of principal concern was the relative merits and demerits of L-band vs. VHF 
frequency. Assistance was also given in explaining differences in statistical processing of 
scintillation data by AFCRL and COMSAT. 

In summary, the Fellowship Program afforded a valuable opportunity for research 
and exchange of ideas and viewpoints. 
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ELECTRON IMPACT IONIZATION OF HYDROGEN NEAR THRESHOLD 


Yiikap Hahn 
Department of Physics 
University of Connecticut 
Storrs, Connecticut 


The study reported here has been carried out in collaboration with Dr. A. Temkin 
and Dr. A. Bhatia of the Lab. for Space Studies, NASA. 

When electrons of about 1 Ryd. energy collide with a hydrogen target, there are an 
infinite number of channels opening up, and reactions of the following forms are possible: 


ej + (e‘ + p+) - e. + (e" + 

- e- + (e- + p+)__ 
^ e j + e" + p'^ 


elastic 


inelas tic 


ionization 


Theoretical analysis of this multichannel problem is complicated further by the long- 
range character of the interactions involved and also by the presence of three-particle 
(ionization) channels. 

The ionization problem near the threshold has been discussed by many people in 
the past, and these studies indicate the correct threshold behavior of the ionization cross 
section depends critically on the dynamical screening of the proton charge by one of the 
electrons as they come away from the proton. The first ingenious classical argument 
of Wannier^ gave the yield Y to depend on E as E * , while the assumption of complete 

screening2 gives the behavior Y(E) ~ E‘ ® as E - 0 and the case of no screening ^ Y(E) ~ 

E More recent investigations ^ of the problem using the WKB-type wave functions 
seem to reproduce the Wannier result. In all cases, however, the most probable con- 
figuration of the two electrons as they move away from the proton is given by 6^^ = n , 
i.e., the electrons are coming out in the opposite directions. 

In an attempt to treat the problem completely quantum mechanically, we have 
introduced a model in which the interelectron potential e^/r ,2 is replaced by e^/(rj + rj), 
which is in fact exactly the same as the original potential for = n. This new poten- 
tial simplifies the analysis a great deal since, e.g., for the zero total partial wave L = 0, 
each-fj =1^ component uncouples. Andyet, the all-important dipole terms are still 
present, as - r'* + (r^ + r^)-! = - v^/t \ for Tj » r^. 

The general formalism adapted for this problem is the projection method of Fesh- 
bach® for the low energy problem at E = -n"^ (Ryd.), where n is the principal quantum 
number and large. One of the motivating reasons for this approach is the earlier calcu- 
lation by Temkin, Bhatia and Sullivan^ in which they showed that the resonance states 
arising from the degeneracy of levels for each n and different 'f ^ n are the dominant 
components of the optical potential and can be represented by unscreened Coulomb wave 
functions. This in principle allows one then to construct a reliable potential. 
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We have analytically solved the problem in the one- level approximation, and found 
that Y(E) ~ E for this case, with or without the one-term optical potential. However, there 
are of course an infinite number of levels contributing to the potential, and thus the result 
has to be generalized to include aH the levels. Their inclusion requires a detailed examin- 
ation of the dependence on n of the widths and shifts of these levels. Preliminary indica- 
tions are that their total contribution may not change much the one- level result insofar as 
the threshold behavior is concerned, but it is not yet clear what the precise magnitude of 
the change would be. We have not completed the study of this model, but we feel that many 
of the conceptual difficulties in the formulation of the theory have been cleared up by the 
present study. We hope to continue the work on the model and also to be able to apply the 
approach to the real physical case. 
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ANALYSIS OF DATA FROM THE FILTER WEDGE 
SPECTROMETER (FWS) ON NUMBUS D 


Edward A. Harms 
Department of Physics 
Fairfield University 
Fairfield, Connecticut 


The summer was spent investigating data received from the Filter Wedge Spectrom- 
eter (FWS) on Nimbus D. The spectrometer scans the wavelength interval from 1.26 to 
2.45 microns. The general idea was to determine how much useful information (from a 
remote sensing point of view) is contained in this interval. 

Two subprograms were investigated. The first was the determination of ground 
truth data available. The second was to determine the information content in this inter- 
val. The quality of the ground truth data determines to a great intent the possibilities 
for success in the second problem. 

Various methods were developed to aid in the interpretation of the data. As a re- 
sult the general qualitative natxire of the spectra are felt to be fairly well understood. 

Some preliminary results were also obtained from a statistical analysis of the spectra. 
Results were encouraging although inconclusive due mainly to the lack of sufficiently 
detailed ground truth data. 
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X-RAY STUDY ON MATERIALS 


Ernest H. Henninger 
Department of Physics 
DePauw University 
Greencastle, Indiana 


Some techniques and theories related to X-Ray studies of materials were studied. 
The work fell in four categories: (1) Preparing and examining powder samples from a 
nickel- titanium alloy; (2) Learning the principles and operating techniques of a Buerger 
precession camera; (3) Studying the techniques of X-Ray stress topography; (4) Studying 
the Kinematical and Dynamical Theories of X-Ray diffraction. 

(1) Powder samples from the bulk ingot were readily identified as being a Nizti 
alloy by indexing the forword reflection Debyescherrer- lines. Of more interest was the 
composition and structure of small fibers which formed on the surface. Attempts to re- 
move individual fibers by etching and cuttii^ failed, but filings from the surface were 
examined. Also back reflection lines from the fibrous surface were recorded. All pat- 
terns showed the Ni^Ti series of lines^‘^ with no evidence of the pure Ni or NiTi phases. 
Only one observed line was unaccounted for. The only metal which could fit both that 
line and a companion group consistent with the other experimental lines is Josephenite, 
a nickel-iron alloy. The sample alloy is NijTi with some possible iron contamination. 

(2) The direct goal was to learn how to take and index precession photographs. 
These are direct magnifications of reciprocal lattice planes.^ Confusing photographs 
early on led to closer examination of the physics of the effects. It was found that usable 
cone- axis and upper level precession photographs could not be obtained using copper 
radiation unless the associated crystal axis has a parameter greater than 10 A°. 
Molybdenum radiation was not available. Hence the early samples were replaced with 

a CaW 04 crystal with large lattice constant, c = 11.4 A°. Although a very thin slice 
was cut the absorption was too high to give good pictures. Nevertheless the symmetry 
of the crystal and the lattice parameters were confirmed. The study did teach the 
method and emphasized the requirements of proper wavelength, crystal size and ab- 
sorption, camera settings, film type and exposure, and crystal alignment. 

(3) X-Ray Stress topography covers a number of techniques^ for detecting 
individual dislocations in perfect or slightly strained crystals. Such stdies have become 
valuable with the availability of perfect crystals for semiconductor devices technology. 
Defects introduced during processing are known to degrade the performance of such 
devices 

Defects are made visible either by the direct image contrast due to kinematical 
scattering from disturbed lattice planes in the vicinity of a defect, or by dynamical ef- 
fects, including the disruption by a defect of anomoulous scattered inten sity'^®\ The 
author hopes to initiate work in this area in connection with his colleague. 

(4) The X-Ray diffraction theories^®^ were studied in order to better understand 
the precession method and the stress topography. The more familiar kinematical theory 
is known to be fundamentally incorrect, but it is useful nonetheless for structure de- 
terminations, since most crystals have a mosaic structure. Scattering from imperfect 
crystals, is best understood in terms of the Ewald construction. In a perfect crystal 
multiple scattering must be taken into account, and the resultant energy transfer pattern 
through the crystal is determined by the complex electromagnetic wave field which can 
exist, subject to Maxwell's equations and the electron distribution. The possible waves 
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are still related to those reciprocal lattice points very near the Ewald sphere, but the 
single Laue point is replaced by a set of points located on sheets of a dispersion surface. 
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A QUASI-OPTIMAL CONTROL SYSTEM 
FOR STEERABLE STRUCTURES 


Eugene W. Henry 

Department of Electrical Engineering 
University of Notre Dame 
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2. INTRODUCTION 

The precise and reliable conti'ol of mechanically steerable antennas and telescopes 
under a variety of environmental conditions is essential for the successful operation of 
the NASA Satellite Tracking and Data Acquisition Network (STADAN) and the Manned 
Space Flight Network (MSFN). STADAN has 73 steerable antennas in 16 locations 
throvighout the world and MSFN has 47 steerable antennas at 21 sites \ 

GSFC directs continuing in-house and contracted research on control systems for 
antennas and telescopes. Recent supported research includes the application of on-line 
optimal digital control^, and model- reference adaptive control'* to the antenna control 
problem. A significant in-house development was a precise digital- analog control sys- 
tem for an optical mounts This moxmt, consisting of a two-axis drive for a 24-inch 
diameter telescope 10 feet long, was available for experimentation with the quasi-optimal 
control system as described in this report. 

This work consisted of three concurrent interrelated efforts; development of a 
model of the plant (system to be controlled) from theory and from experimental data; simu- 
lation of the model and the quasi-optimal control on an analog computer at NTTF; and im- 
plementation of the control system on the optical mount at GORF. 


3. STATEMENT OF THE PROBLEM 

The problem is stated in terms of controlling one axis of the 24-inch telescope at 
GORF, but similar specifications obtain for steerable structures in general. The slew 
mode is used for acquisition of a target, after which the track mode is employed for pre- 
cise following. The controlled system must meet the following specifications: 

Maximum acceleration: 2 deg/sec ^ slew, 0.2 deg/sec^ track 

Velocity range: 5 deg/sec maximum slew, .0005 to 

5 deg/ sec track 

Tracking accuracy: ±1 arc second 

These requirements are to be maintained with noise of specified properties on the 
tachometer and error signals, and with variations of X-axis inertia between 1900 and 
4600 Ib-ft-sec^ caused by motion of the Y-axis. 
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4. THE QUASI-OPTIMAL CONTROL 

In the slew mode, it is desired to acquire the target in minimum time, subject to 
the constraints on maximum acceleration and velocity. It is well known that the time 
optimal control of a double-integral plant is accomplished by means of a relay servo 
with a quadratic switching curve. Time optimal control of higher order plants, although 
theoretically possible, require extensive on-line computation. 

The quasi- optimal control method employs appropriate compensation networks be- 
tween the tachometer and the servo amplifier so that the resulting minor loop has a trans- 
fer function approximating a double integration. The outer loop contains a relay controller 
with the quadratic switching criterion. Velocity limiting is provided by a dead-zone feed- 
back circuit from the tachometer in the minor loop. 

In response to a step input in the slew mode, this sytem will cause the mount to 
accelerate at 2 deg/sec^ until the velocity reaches 5 deg/sec, and maintain that velocity 
until the controller commands an acceleration of -2 deg/sec^ to bring the mount to the 
commanded position in minimum time without overshoot. 

In the track mode three gains are changed simultaneously in the controller to reduce 
the maximum commanded acceleration to 0.2 deg/sec^ and to change the dynamic operat- 
ing region of the quadratic switching curve circuit. 

Advantages of this quasi-optimal controller are. 

a. It is relatively simple to implement with either analog or digital devices. 

b. It essentially provides time-optimal control for step inputs in the slew mode. 

c. The system will reject torque disturbances within the capacity of the motor. 

d. It has low sensitivity to changes in load inertia. 

e. No transient disturbance is caused by switching between track and slew modes. 


5. DEVELOPMENT OF THE MODEL 

The model was developed from theory and experimental data as a fifth order plant 
to include an oscillatory structural mode as well as the electrical and mechanical 
aperiodic modes. Static friction, coulomb friction, and power amplifier saturation were 
also included in the model. Methods were developed to readily determine the values of 
friction, inertia, and other plant parameters from simple experiments on the mount. 


6. SIMULATION AT NTTF 

The model of the plant and the quasi-optimal controller were simulated on an 
Electronics Associates, Inc. TR-48 Analog Computer at NTTF. The simulation included 
provision for studying the effects of random disturbances on the tachometer signal, the 
error signal and the output (torque disturbance). This simulation was used to determine 
optimal settii^s of controller parameters and to investigate the behavior of the system 
in both track and slew modes with various disturbances and reference inputs. 
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7. IMPLEMENTATION AT GORE 

The controller circuitry was implemented on another EAI TR-48 Analog Computer 
at GORE. The tachometer and position-potentiometer signals from the telescope were 
connected to the computer, and the computer provided the driving signal to the servo 
amplifier. 


8. RESULTS 

The results of the computer simulation indicate that the quasi-optimal control is 
feasible and can meet the specifications while providing the advantages listed in Para- 
graph 4. Limited experiments at GrORE verified the desired system behavior in the 
slew mode and the absence of transients in switching between track and slew modes. 
Performance of the telescope in the track mode was not satisfactory due to excessive 
noise. 


9. SUGGESTIONS EOR EURTHER RESEARCH 

a. The sources of noise in the physical system should be isolated and analyzed so 
that appropriate filters can be incorporated. 

b. Implementation of quasi-optimal control by digital rather than analog techniques 
should be investigated as a means for providing precise control while eliminat- 
ing much of the drift and noise which plagues the analog controller. 

c. Alternative approaches to the design of compensation for the inner loop and 
the possible inclusion of another minor loop enclosing the relay element to 
provide a system which is stable for all values of gain^ should be investigated. 
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A FEEDBACK COMMUNICATION SYSTEM 
USING CONVOLUTIONAL CODES 


Joseph L. Katz 

Department of Electrical Engineering 
Purdue University 
Lafayette, Indiana 


In recent years the application of feedback to communication systems has received 
a generous amount of attention in the literature and certain theoretical aspects of feed- 
back systems has been studied in detail. Two-way communication systems have the 
capability of transmitting information about the current state of a message being decoded 
at the receiver back to the transmitting point. The returned information can be used to 
simplify the coding and decodii^ operations in the forward channel and to provide a lower 
probability of error for a given coding delay than could be achieved without feedback. A 
potentially useful application of feedback communication systems is in the design of ef- 
ficient data retrieval systems for space vehicles, where the transmitter power of the 
spacecraft is restricted to be several orders of magnitude less than the transmitting 
power of the based equipment. 


In general communication feedback systems can be classified according to the type 
of information carried by the feedback channel. The term "decision feedback" is used 
when the feedback channel is used only to report the decision of the receiver as to the 
acceptability of each received message unit, and "information feedback" when the feed- 
back channel is employed to report information about each received message unit to the 
transmitter, the decision to accept, or reject and correct, subsequently being made at 
the transmitter. 

In either type of system once a decision is made that the message cannot be in- 
terpreted, the transmitter may repeat, in whole or in part, the ambiguous message unit 
either usii^ the same code and transmitter power or by agreement, shifting to a higher 
redundancy code and transmitter power or any combin^ion of the two. The adaptive 
nature of feedback systems permit the efficient matching of the signals to variations in 
the channel conditions. This in turn permits the setting of an upper limit on the proba- 
bility of error, whereas for a imidirectional link the signals must be chosen such that 
the error probability will not be excessive at the poorest signal-to-noise ratio which is 
anticipated. 

During the fellowship period a decision feedback communication system, using 
convolutional encoding in the forward channel, was studied. Results were obtained for 
the case of sequential decoding and partial results for maximum likelihood decoding of 
a 1/2 rate systematic convolutional code. 

The feedback system incorporating convolutional encoding- sequential decoding used 
a simple request repeat strategy in which repeats are requested for frames in which a 
computation overflow occurred. The repeated frames are retransmitted using a 1/2 
rate non- systematic code and a 1/4 rate code is formed from the first and retransmitted 
information. Decoding is attempted a second time using the resulting 1/ 4 rate code. It 
an overflow occurs the second time a detected error is declared. A full scale computer 
simulation of this sytem revealed substantial improvement in the computational overflow 
rate over the unidirectional system for E^/No < 2.0 db. 
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The problem of computational overflow does not exist for maximum likelihood de- 
coders but decoders of this type can only be built for short-constraint ler^h convolu- 
tional codes. This limits the amount of improvement in bit error over no coding which 
can be achieved. A feedback system similar to the one described earlier was imple- 
mented using a maximum likelihood decoder. Repeats are requested for frames in which 
the metric falls below a given threshold. Indications are that this system will improve 
the bit error rate of the short constraint length code; however, results are sketchy at 
this time. Work on this topic will continue under a NASA grant. 
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INVESTIGATION OF WALSH FUNCTIONS FOR 
SIGNAL PARAMETER EXTRACTION 


Thomas F. Krile 

Department of Electrical Engineering 
Rose Polytechnic Institute 
Terre Haute, Indiana 


I. INTRODUCTION 

As signal processing systems consume ever increasing amounts of data, much of it 
having low information content, the problems of bandwidth reduction and signal parameter 
extraction become more important. Spacecraft in particular are required to gather large 
amounts of data and transmit it to earth. Due to the inherent power restrictions aboard 
spacecraft, and the complex coding required to minimize errors, it becomes apparent that 
some signal processing should be done on board so that only the most significant signal 
parameters are transmitted. 

Other signal processing systems, such as in ECM or surveillance, require essen- 
tially real time data processing. In such cases, it is useful to examine only a few speci- 
fic signal parameters, and the problem is to determine which parameters are the most 
significant and most quickly and easily obtained. 

One well-known way to analyze signal parameters is to use a given set of basis 
functions to transform the signal to another representation. Then the new signal space 
is examined to see if any important information is more apparent than it was in the 
original space. A typical example of this process is the Fourier transform analysis of 
signals. The resulting frequency spectrum has long been used in single-dimensional 
signal analysis, and more recently has shown itself to be a powerful tool for two- 
dimensional image processing^ (i.e., parameter extraction for pattern recognition). 

The Fourier transform has also been proposed for bandwidth reduction schemes in 
which only a selected portion of the spectrum is transmitted and the signal is recon- 
structed for this truncated spectrvim.^ 

Progress in integrated circuit technology has made the hitherto academic Walsh 
transform into a practical tool for digital signal processing.® With suitable definitions 
for the Dirac delta function and the convolution operation, one can show many properties 
of the Walsh transform which are analogous to those of the Fourier transform. Also, 
relationships analogous to Parseval's theorem and the Wiener- Khinchin theorem can be 
shown to hold. This means that one can perform data reduction and filtering in the Walsh 
domain. 

The first part of this report presents definitions and properties of the one- and 
two-dimensional Walsh functions and Walsh transforms. The primary part of this paper 
is devoted to uses of the Walsh spectrum of two-dimensional images for parameter 
extraction. A set of "correlation masks" is shown which, when applied to the image, 
result in the individual elements of the Walsh transform. The same masks can be ap- 
plied to the transform to recover the original image. 

As in the Fourier transform case, the orientation and separation of line structure 
in an image can be determined from the Walsh transform. However, there is an angle 
ambiguity, since all the amplitude information in the Walsh plane is contained in one 
quadrant, whereas in the Fourier plane it is in two quadrants. 


69 



The most troublesome aspect of the Walsh transform is that it is not translation 
invariant, as the Fourier transform is. A form of power spectrum is developed which 
is translation invariant and which reduces the number of spectral data points from 
(2’’)^ to (p + 1)^ where 2" is the order of the Walsh matrix. This power spectrum pre- 
serves the angle information contained in the original transform. Its other pattern 
recognition aspects have yet to be investigated. 


II. DEFINITIONS 

Walsh functions form a closed set of orthogonal functions, bounded in time (or 
space), which can be used as bases for signal representation.® In this report, the dis- 
crete version of the Walsh functions is used. 

Each Walsh function can be viewed as a row vector Wal(i) of dimensionality 
N = 2” where p is an integer. Let Wal(i,m) refer to the m*’’ component of the i‘'’ 
function, i = 0, N - 1 and m = 0, N - 1. The vector components consist of +l's and 
-I's, so the functions are binary in nature, and Wal(l), is (1,1,1, . . . , 1) by definition. 

Walsh functions are ordered according to the number of zero crossings (transi- 
tions from ±1 to +1) in the interval of definition. This gives rise to the generalized 
frequency concept of "sequency" where, if j > i, then Wal(j), will have a greater 
sequency (greater number of zero crossings per unit interval) than Wal(i). Using 
this ordering, the Walsh transform sequency plane will be roughly analogous to the 
Fourier transform frequency plane. 

The functions themselves may be defined via several routes. They may be written 
as products of Rademacher functions taken in a certain order.’ They may also be ob- 
tained as the solutions of a difference equation'* (as the Fourier sine and cosine com- 
ponents arise from a differential equation). We consider here a function generating 
method proposed by Swick** which allows one to obtain individual functions independently 
by taking advantage of the sequency ordering and the symmetries inherent in the func- 
tions. 


The first step in forming Wal(i) is to express the sequency number i in binary 
form as b^bj . . . bp_j where the decimal value of i is given by 


i = bg 2P-* + bj 2P'^ + . . . + bp., 2®. 


( 1 ) 


We let Wal(i,0) be +1 for all i. Now if b^ is 1, the m = 0 and m = 1 components of Wal(i) 
are skew symmetric, otherwise they are symmetric. If bj is 1, the (0,1) and (2,3) com- 
ponents of Wal(i) are skew symmetric, otherwise they are symmetric. If b 2 is 1, the 
(0, 1, 2, 3) and (4, 5, 6, 7) components of Wal(i) are skew symmetric, otherwise they are 
symmetric. In general, if b^ is 1, the (0, 1, . . . , 2*‘ - 1) and (2’', . . . , 2''*^* - 1) com- 
ponents of Wal(i) are skew symmetric, otherwise they are symmetric. An APL function 
using this generating method is included in the Appendix. 

As an example of this procedure, consider a system of order N = 8 (p = 3) and the 
function Wal (6). First we have i = 6 = 110 = bgbjbj. Since b^ = 1, we find 1 = Wal (6, 0) 

= - Wal (6, 1). Also, bi = 1, so Wal (6,0) = -Wal (6,3) = 1 and Wal (6,1) = -Wal (6,2) = -1. 
Finally, b^ = 0, so Wal (6,0) = Wal (6,7) = 1, Wal (6,1) = Wal (6,6) = -1, Wal (6,2) = Wal 
(6,5) = 1, and Wal (6,3) = Wal (6,4 = - 1. Thus we have Wal (6,) - (1, -1, 1, -1, -1, 1, -1, 1). 
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An important propGrty of the Walsh functions is that they are orthogonal, i.e. 


N-I 


E 

isn 


Wal (i, j) Wal (k, i) = 



i ?^k 
i = 1. 


(2) 


In order to obtain the Walsh transform, a matrix is set up composed of rows which 
are the individual Walsh functions in ascending sequency. As an example, consider the 
case for N = 2® = 8. The Walsh matrix is: 

Sequency 

“l 1 1 1 1 1 1 1~| 0 

1 1 1 l-l-l-l-l 1 

1 l-l-l-l-l 1 1 2 


W = 




1 

1 


1 


1 

- 1 
- 1 


- 1 - 1 
- 1 1 

- 1 1 


1 1 - 1-1 

1 - 1-1 1 
_1 1 1-1 


1-1 1 - 1-1 1-1 1 


3 

4 

5 

6 


(3) 


1 - 1 


- 1 1-1 1 1 


7 


Note that the matrix is symmetric, it is in fact a re-ordering of a symmetric Hadamard 
matrix (a related signal processing operator is the Hadamard transform*^ ), 


If the signal to be transformed is an N-dimensional column vector X = (Xj, x^, 
. , , . X ), then the Walsh transform is defined as 

" N ' ” 


f = (1/N))K‘X=(1/N)){[X 


(4) 


where t denotes the transpose. The inverse operation is given by 


X = W T. 


(5) 


When the signal is an N‘*’ order square matrix X, the transform is found from: 


T = (1/N1 W‘XW-(1/N))KXW. 


( 6 ) 


The inverse is then 


X = (1/N) WTW, 

and we see that the two-dimensional Walsh transform is self-inverse, 
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(7) 



m. PROPERTIES 


By maMng suitable definitions, one can show several properties of the Walsh trans- 
form which are analogous to those of the Fourier transform,® These properties can give 
us a hint as to what features or parameters will be significant in the sequency plane. 

The discussion to follow will consider one-dimensional signals and transforms; ex- 
tension to two dimensions is straight-forward. Let a delta function be defined as a vec- 
tor s = (1,0, . . . , 0). Then the Walsh transform of 8 is the vector T = (1/N,1/N, . . . , 
1/N), where N is the order of the system. Thus we see that a delta function in time (or 
space) transforms into a constant in the sequency domain, and vice versa. In general, a 
compression in the time domain results in expansion in the sequency domain. This same 
result occurs in the Fourier case. It should be pointed out, however, that the time- 
sequency bandwidth product is always finite for the Walsh system (because both the sig- 
nal and basis functions must be finite), but the time-frequency bandwidth product is in- 
finite for the Fourier transform. 

Using a special relationship due to Gibbs^, we can define "logical convolution" of 
two sequences (vectors) X and Y as 


N~1 

\ (s) = g/N)^ X(i) Y(s 

i-0 


i); s = 0, N - 1; 


( 8 ) 


where ® indicates addition modulo 2, We abbreviate this to = X © Y. Now letting the 
Walsh transforms of X and Y be > respectively, and using the addition formula 

for Walsh transforms (Wal(i,m) Wal(j,m) = Wal(i © j, m)), we have 


N*1 

(s) (i) Ty (i) Wal (i, s). 

iso 


( 9 ) 


Thus we obtain the logical convolution of two functions, X and by takii^ the_Walsh 
transform of the term-by-term (scalar) product of their transforms, T^ and T^. The 
same procedure is applied to find the normal convolution of two functions using Fourier- 
transforms. 

Using the definition of logical convolution, we can define the logical autocorrelation 
of a function X as 


= X © X. 


(10) 


Also, we can define a Walsh power spectrum for X where 


s,. (i) [T^ i = 0, N - 1. 


( 11 ) 


This is analogous to the Fourier power spectrum definition. Now usiTO Eq, (10), we find 
that the logical autocorrection of X and the Walsh power spectrum of X are a Walsh 
transform pair. This statement is the analog of the Wiener-Khinchin theorem in the 
Fourier case. 
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The results given above can also be combined to show a Walsh analog to Parseval's 
theorem, i.e. 


N~ 1 N~ 1 

a/N)^^ [X(i)]^ = Y' [T^ 


i=0 


j*0 


())]' 


( 12 ) 


These Walsh transform analogs of the familiar Fourier transform theorems justify the 
investigation of the Walsh spectrum for filtering, bandwidth reduction, and parameter 
extraction operations. 


IV. PARAMETER EXTRACTION 

A set of correlation mask operators has been found which is useful for determining 
parameter extraction properties in the Walsh sequency domain. Consider a matrix mask; 

, consisting of +l's and -I’s, which is of order N. The mask is obtained by modulat- 
ing the column vector Wal (i) by the row vector Wal(j). An example of the masks , 

, . . . , ]^33 for a fourth order Walsh system is shown in Figure 1. 

Let Q = M X 5^ be the scalar (term-by-term) matrix product of mask M , and the 
image matrix X. We find that ^ ‘ ' 


N«1 N - t 

(1 y~'c (n. m) = T (i, j) (13) 


where J = (1/N) ^ Thus, using 11(1^. as a template on the image, and summing the 

elements of the resultant matrix scalar product, we obtain the (ij) term of the image's 
Walsh transform. If the image pattern correlates highly, either positively or negatively, 
with the mask M, . pattern, the (ij) term of the transform will be high. If the correlation 
is slight, the corresponding transform term will be small. This means that the patterns 
embodied in the masks can be studied to determine what image parameters will be de- 
tectable in the sequency plane. 

Since consists of all +l's, we observe that the transform's T(0,0) terms must 
contain the average, or DC, value of the image X. This parameter is useful both in pat- 
tern recognition and statistical signal processing. 

A study of the masks also tells us how lines or edges in the image plane appear in 
the sequency plane. A vertical line in the image plane will correlate most highly with 
masks , M 02 . etc., and thus will appear as a horizontal line in the sequency 

plane. Similarly, a horizontal line in the image plane transforms to a vertical line. In 
general, a line at ar^le 6^ {e^ in the image plane will transform to a spectrum whose 
largest magnitudes fall along a line of angle -4^ in the sequency plane (see Figure 2) with 
the following relationships: 


for 0 < < - 90'^: ^ - (90° + 9 ), 

(14) 

for - 90° < < _ 180°; 1 /; = 90° + 
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Figure 1 , 


This gives rise to an angle ambiguity. If ^ is known, the image plane angle ' is either 
-(>/' + 90°) or ^ - 90°, i.e. 0 is a double-valued function of v. A slight rotation of the 
image could be used to resolve the ambiguity in practice. 

It should be noted here that the correlation masks can be applied to the transform 
plane and used to obtain the individual elements of the image, since the transform has 
been shown to be a self-inverse operation. Thus any algorithm or physical device to 
use these masks will be bi-directional. 
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T(O.O) 


Image Plane X Sequency Plane T 

Figure 2. 

One of the primary operations done in the frequency plane is filtering, and we 
would like to extend the concept to the sequency domain. As a common example, con- 
sider the application of high-pass filters in the frequency domain to bring out lines and 
ec%es in the reconstructed signal. The analogy between sequency and frequency leads 
to the idea of trying the same thing in the sequency domain. Here we run into a primary 
difference between the Walsh and Fourier transforms. High pass sequency filtering will 
indeed bring out edges, but the edges which are enhanced depend on where the image is 
located in the im^e plane. Thus, while the Fourier transform is translation invariant, 
the Walsh transform is not. If the image is moved vertically or horizontally, all the 
values in the sequency plane will change, in general. 

In an attempt to get around this problem, a two-dimensional power spectrum was 
evolved which is translation invariant and which reduces the number of spectral data 
points from (2^)^ to (p + 1)’ where N = 2’’ is the order of the Walsh matrix. The power 
spectrum is obtained via a matrix operator which sums the energies of certain frequency 
groups, namely the odd harmonics of 1, 2, 4, 8, 16 etc. cycles (not zero crossings). If 
we let y = X X gr be the scalar matrix product of the Walsh transform matrix, with itself, 
the power spectrum P can then be defined as: 




p = A v A‘. 


(14) 


Here' the matrix operator 4 is of order p + 1 by N. As an example, for N = 16, the matrix 
^ is defined in the function INPWR16 in the Appendix. 

This power spectrum is translation invariant and preserves the average value and 
line angle information of the original transform. Its other pattern recognition aspects 
have not yet been studied. 


V. CONCLUSIONS 

In this paper, several properties of the Walsh transform have been noted which, by 
analogy to Fourier transform properties, show that the Walsh sequency plane has poten- 
tial uses in filtering and bandwidth reduction schemes. 
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The pattern recognition aspects of the sequency plane were also investigated. It was 
shown that lines and edges in the image plane can be detected in the sequency plane. A set 
of "correlation masks" was developed which one can use to find certain other image pat- 
terns which can be identified in the Walsh spectrum. Since the Walsh transform is not 
translation invariant, it is difficult to get results which hold for arbitrary patterns. A 
type of power spectrum was evolved which is translation invariant and which preserves 
some of the parameter extraction properties of the original transform. 

The characteristics reported here, along with the features of simplicity and efficient 
representation of digital signals, show the Walsh transform to be an important candidate 
for signal processir^ applications. A technique has recently been proposed ‘ for obtainii^ 
the Walsh transform of an image optically, and as other physical implementations for tak- 
ii^ the transform are devised, its use in signal processing is certain to grow. 
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APPENDIX 


This portion of the paper contains the main APL programs which were used during 
the study. Some of them have been referenced in the text as examples of the operations 
performed. The whole set of functions is included to provide continuity and to show the 
simplicity and utility of the APL system for vector and matrix calculations. The first 
functions were written and running after only two day's introduction to the APL language 
and terminal. The functions are not optimized in terms of length or time, but they get 
the job done with a minimum of experience required. 

INIT16 initializes the variables needed to work with Walsh functions of order 16. 

N = 2’’ is the order, M = p, and V is a row vector of length M, composed of 2's. 

WALSH is the function which generates the Walsh matrix, IWALSH, of order N, 
according to the algorithm given in the text. 

XFRM takes the Walsh transform of a square image matrix P of order N. The 
transform appears as the matrix TRANS. 

FILTER initiates a filter mask matrix composed of all +l's. Appropriate elements 
are then set to zero to provide the desired filterii^. 

NEGINV takes the Walsh transform of the filtered sequency spectrum and recon- 
structs the filtered image amplitude. 

INPWR16 initiates the matrix operator for obtaining the power spectrum (as 
developed in part IV) of an order 16 system. POWER16 then finds the power spectrum 
and displays it. 
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APPENDIX 


APL Functions For Walsh Image Processing 
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PROPAGATION DELAY IN THE ATMOSPHERE 
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ABSTRACT* 


The structure of the lower atmosphere and the ionosphere 
are discussed and an introduction to ray tracing the atmosphere is 
given. Models for each region appropriate for satellite-to-earth 
communication in the frequency range 1 GHz are discussed. As- 
suming a spherically symmetric atmosphere, propagation delay is 
computed for the lower atmosphere and for the ionosphere. The 
total atmospheric delay is also computed as a function of atmospheric 
conditions, frequency and altitude of the signal source. 
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*This work was published in X-521 -70-404. 
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CURRENT INSTABILITIES IN OXYGEN-DOPED 
GALLIUM ARSENIDE 


Alfred George Lieberman 
Department of Electrical Engineering 
University of Maryland 
College Park, Maryland 


RESEARCH SUMMARY 

Electrical current passing through prepared gallium arsenide was observed to 
display large amplitude, low frequency oscillations when suitably illuminated and 
properly refrigerated. Oxygen-doped gallium arsenide samples having a room- 
temperature, dark resistivity of 8.6 x 10“' n cm were found, at 100°K, to exhibit at 
least two distinct modes of oscillation. The threshold electric field required to sus- 
tain oscillation is determined by the wavelength of the incident light and the mode of 
oscillations. The lower threshold mode exhibits a sinusoidal current waveform. At 
larger electric field intensities a second threshold exists at which the current wave- 
form abruptly changes into a repetitive spike of slower periodicity. Pronounced photo- 
excitation is observed for wavelergths between 0.9 and 1.6^i, corresponding to the 
energies separating the conduction band from the valence band and the oxygen dopant 
levels. For a given mode, the threshold field is observed to depend upon the wavelength 
of the incident light; the frequency of oscillation is determined by the intensity of the 
light and the length of the crystal sample. Threshold field inten sities vary from a few 
tens of volts/cm for the sinusoidal mode to several hundreds of volts/cm for the spike 
mode. Crystal samples were cut and polished to dimensions typically 0.5 mm x 1 mm 
X 2 mm. 
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SPACECRAFT POWER SYSTEMS 


Frederick H. Morse 
Department of Mechanical Engineering 
University of Maryland 
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My research colleague at NASA was Mr. Joseph Epstein, Head Advanced Power 
Sources Section, Code 761. Two projects were selected for study. 

One project was to develop a method to predict the change in the specific power 
of the various satellite energy conversion systems, presently being used or under 
various stages of development, as the satellite power level increases. In this study 
specific power is defined to be the orbit average regulated system weight. A set of 
relationships were to be developed that would show, for any energy conversion system, 
how the weight of the power system changed as the load requirement increased. The 
predicted variation in specific power with power level was to be compared with data 
from actual satellite power systems. A survey of satellite power systems was con- 
ducted to provide such data. 

The analysis, covering solar cell, isotope, thermoelectric, thermionic, magneto- 
hydrodynamic and dynamic cycle energy conversion systems, has not been completed. 
The survey of satellite power systems, all utilizing solar cells as the primary energy 
conversion device, has been completed. Of the approximately one hundred twenty non- 
military satellites launched since 1959 sufficient power system data was found for only’ 
fifteen. Additional data on five satellites presently being built were also obtained. 
Several preliminary observations can be made. 

1. The maximum specific power is 1.5 watts/1 bm. 

2. For earth orbiting satellites (with the exception of the Intelsat satellites) a 
maximum specific power of approximately 1.1 watts/1 bm occurs at the 100 
watt power level. Low power satellites (20-40 watts) have a typical specific 
power of 0.6 watts/1 bm while high power satellites (100-300 watts) have a 
typical specific power of 0.9 watts/1 bm. 

3. Over the last ten years the specific power has gradially increased. Satellites 
in any series, such as the Nimbus series, always show an increase in specific 
power with power level. 

The second project was to investigate the feasibility of a hybrid thermoelectric- 
thermionic power system. In order to prevent back emission, resultii^ in a decrease 
in output power, the anode of a thermionic generator must be cooled. At present this 
cooling is accomplished by heat rejection to space. In the hybrid system this heat will 
pass thru a thermoelectric generator thereby converting a portion of this rejected 
energy into electrical energy. The analysis of this system is continuing. 
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ELECTRONS IN AURORA 
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Victor J, Newton 
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and 
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An experimental program for measurii^ auroral atmospheric emissions in the in- 
frared, visible, and ultraviolet wavelength regions has been developed by D. F. Heath at 
this Center. The program includes ground based as well as rocket and satellite instru- 
ments, Of particular interest is Heath’s 1 meter Ebert-Fastie spectrometer, which is 
capable of rapid scanning at high resolution over an auroral band, as viewed from the 
ground. This is a rather useful measurement to make since it has been shown that the 
temperature of the emittii^ region can be reasonably well inferred from the data. For 
this purpose the N^ 3914 A (0,0) band is especially weU suited. 

A complete theoretical analysis of auroral phenomena involves detailed knowledge 
of the density and temperature structure of the atmosphere, the electron- and to some 
extent ion- impact cross sections for exciting each state, the magnetic and electic field 
components, and finally the incident particle flux spectrum. For most auroras, electrons 
are the dominant source of excitation as they degrade from their initial incident energies. 
Thus, for example, in order to predict the height profile of the 3914 A (0,0) band, one 
must know how electrons deposit their energy as a function of altitude. 

The initial work in the overall problem has been to find a simple way of correlat- 
ing a given incident electron spectrum to the resulting 3914A emissions and to examine 
how these emissions, integrated along a line of sight are seen from the ground. For 
simplicity, initially we have assumed that a mono- directional steady-state beam of elec- 
trons is incident on a layer of gas, thin with respect to the total range of the electrons in 
the gas. After passage through the gas, the electrons will have lost some of their energy 
through excitation and ionization processes, and what emerges is a flux distribution. For 
the most part, inelastic collisions are predominantly forward scattering; thus, we assume 
no deviation from a straight line path. Although this is not completely realistic it is 
hoped that phenomenological corrections can be brought in to produce agreement with 
experiment and with more detailed Monte Carlo calculations. As a test, the results 
should be somewhat consistent with a Landau distribution. 

We consider the energy spectrum to be divided into a series of bins of specified 
width AE, and write the density Nj(r) of electrons in the ith bin, as a function of the 
position r. We expect this density to follow a continuity equation including a source term 
to account for electrons being deposited in the ith energy bin through deposition from 
higher energy bins, and a loss term to account for deposition from the iUi bin to lower 
energy bins. We can write, therefore: 


3 N. (7) 

3 r 


= V. VN. (r) = 


N. (7) 

L (r) 


+ 


Z— I T: (t) 

j<i ’ 


A cr. . 

j_± 


a. 

j 


( 1 ) 
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where we have labeled the bins consecutively such that j = 1 corresponds to the highest 
energy bin. The collision time (^) of an electron of energy E is related to the mean 
free path = v. t. . In this equation the ratio of partial to total cross sections A cr. ./o-. 
represents the probability that an electron will go from bin j to bin i. ‘ 

Considering a one- dimensional steady state solution for a mono- energetic beam of 
electrons impinging on a slab of gas of density n(x), the continuity equation simplifies to 


d N. (x) N. (x) ^ — , Vj N, (x) ^ CTj i 

d X \ (X) \ (X) O'. 

i<i ‘ J J 


where represents the mean free path for particles in bin i. Since (x) = 1/a, n(x) 
and since the electron flux Ij (x) is just Nj (xjv^ , it is convenient to translate Eq. 2 'into 
the flux equation; 


d r (y) 

d y 


- I (y) cr. 


+ ^ I. (y) A a.. 
j<i 


(3) 


in which a change of variable dy = n(x)dx has been introduced. The procedure is to solve 
Eq. 3 to obtain the electron flux, and thereby the volume excitation and ionization rates 
as a function of altitude and incident energy flux. 

Equation 3 can be solved analytically, in the case where no two cross-sections a. 
are exactly equal, by a finite series solution: ’ 


li (y) = 



(4) 


where the coefficient Aj^ satisfy the recursion relation: 


A. . 


I J 



i<k^i 


and the boundary condition: 


(5) 



= li (O') 


(6) 


The advantage of this solution is its simplicity; numerically however, sizeable 
round-off errors can occur which makes it difficult to work with. An alternative solution 
may be obtained by solvii^ for the Green's Function of Eq. 3. Here the solution to the 
source equation 
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( 7 ) 


-^+ a. ) g. (y, y') = S (y - y') 
, d y V 


Subject to the boundary condition of particles traveling in the positive direction is: 


-o-; (y-y ) , 

e , y > y 


g: (y. y‘) = < 


0 . 


y < y 


( 8 ) 


which yields a solution to Eq. 3 of the form 


(y) = 


I; (O') 


j<i 


I’ 


d y' i;. (y') e 


A a. 




( 9 ) 


which can be solved numerically. The advantages of this solution are that the terms 
within the sum in Eq. 9 are all positive, resulting in a negligible round-off error, and 
that the terms in the solution can be interpreted physically. 

The analytical form of the solution to Eq. 3 has been programmed using recent 
cross-sections in phenbmenlogical form. These studies have not yet produced concrete 
results, but preliminary indications are encouraging for many conditions. 
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CODING AND THE SPACEBORNE ATOMIC CLOCK* 


Jacob C. Sabto 

Department of Electronic Engineering 
California State Polytechnic College 
San Luis Obispo, California 


ABSTRACT** 

A satellite-borne atomic clock would be the crucial com- 
ponent in a communications network which could be synchronized 
to an accuracy approaching one microsecond. With such synchro- 
nization accuracy, the technology of time and frequency measure- 
ment would take a leap which, when translated into practical achieve- 
ments, would open new vistas to the solution of problems in such 
fields as submicrosecond world-wide time synchronization, space 
communication, navigation, air traffic control, etc . 

Mr. A. R. Chi has proposed the orbiting of such a clock in a 
report* submitted to NASA Headquarters. We describe in this docu- 
ment a typical implementation of Mr. Chi’s proposal and consider 
one particular aspect of the system, namely coding techniques for 
transmitting information. Before discussing these techniques, it is 
useful to briefly describe the overall system in order to better ap- 
preciate the role of coding within the system. 


blank 




*A proposal for the Support of the Development of a Space Borne Atomic Clock," May 1970, sub- 
mitted to the Advanced Applications Flight Experiments Program Office, NASA Headquarters, 
Washington, D.C. 

**This is being published under X-521 -71 -81 . 
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A RANDOM DISTRIBUTION OF RADIAL MOTIONS 


William R. Thickstun 
Department of Mathematics 
Clarkson College of Technology 
Potsdam, New York 


ABSTRACT 

Two quite unrelated problems motivated the efforts described 
in this report. One pertains to the question of whether a singularity 
in space time can develop from an initially non- singular distribution 
of matter. The other problem pertains to the initial distribution of 
comet orbits about the sun from which certain "capture" theories 
attempt to predict the observed number of short period comets. For 
both of these problems it is instructive to start with a model con- 
sisting of a random distribution of particles constrained, however, 
to strictly radial motion in a central gravitational field. 


PE1CES>ING PAGE BLANK NOT fILMEQ 


*This paper was published under X-643-70-303. 
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ELECTRON PARAMAGNETIC RESONANCE OF AND 

OPTICAL SPECTRA OF IN ZOISITE Caj AI 3 Sij 0,2 (OH) 


Tung Tsang 
Department of Physics 
Howard University 

PRECEDING PAGE BLANK NOT FILMED D. c. 

and 

Subrata Ghose 
Planetology Branch 
Goddard Space Flight Center 
Greenbelt, Maryland 


ABSTRACT* 

Electron paramagnetic resonance and optical absorption spec- 
tra have been used to clarify the local environments of transition 
metal ions and their distributions among various Al- and Ca-sites 
in a gem quality zoisite crystal from Tanzania. The EPR spectra 
due to Mn^'^, Fe^'^ and two t3qDes of have been interpreted by 
the spin Hamiltonian 


;5H-g-S + I-A'S+D 




+ E (S2 - Sp 


in the principal axes system x, y, z. The Mn^'^ ions occupy one of 
the Ca-sites, probably Ca(l), with point group symmetry m; g = 2.003 
± 0.005 and |A| = (85 ± 2) xlO"'* cm"i, bothisotropic and also |D| 

= (103 ± 5) X 10''' cm*', |E I = (34 ± 2) x lO""* cm”', in general 
agreement with other oxides and hydrates . Highly anisotropic hyper- 
fine and small zero-field splittings (| d| < 0.014 cm"* ) are ob- 
served for V^'^, which is the opposite of the usual situation. One 
of the two types of apparently occupies the same site as Mn^'^; 
the other V^'*' ion is in a general position, and probably occupies a 
double minima potential around the Ca(2) site. The Fe®^ ions oc- 
cupy the Aljj site, since the point group symmetry is m; the prin- 
cipal axes, however, are displaced approximately 45° from the NMR 
axes of ^Al. Apparently, the local environment changes when AP *' 
ions are replaced by Fe*'''. For Fe**, we found |d | = 0.14 ± 0.03 
cm ' and |e| ~ 0.1 |d 1. The optical absorption ofV^^ and the 
trichroism of the zoisite single crystal are consistent with the two 
types of Al^^ site symmetries, with crystal field parameters Dq = 
1850 and 1400 cm”' at Alj and Aljj sites. 


*This X644-70-319 to be published in: The J. Chem. Phys. (Spring 1971). 



REDUCTION OF OAO-II STAR 
TRACKER FLIGHT DATA 


Robert E. Wilson 
Department of Astronomy 
University of Sovtth Florida 
Tampa, Florida 


With the help of programmer Sy Lee, I have been attempting to analyze data from 
OAO history tapes to ascertain the dependence of measured star brightness on electron 
temperature telescope temperatures, catalog magnitude, and other parameters. Diffi- 
culties in data retrieval have prevented completion of this project, so we plan to com- 
municate by telephone for the next several months until it’s satisfactorily concluded. 
When it is in the production stage, it will be possible to quickly determine the sensitivi- 
ties of the trackers as functions of any relevant variable. 

I have also been makii^ use of the 360-01 computer for a project which is too 
large for the 360-65 in Tampa, This involves fitting a new binary star model to ob- 
servations of MR Cygni, a well observed eclipsing system. This was done in two 
stages - by intuitive adjustments and by differential corrections, and the final corrected 
set of parameters gives a very good representation of this data. I have been invited to 
speak about this work at the XIV General Assembly of the International Astronomical 
Union, in Brighton, Finland. The paper will be given in late August. > 


^preceding page 
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DIGITAL SIMULATION OF FCFB LOOP OPERATION 
IN RFI AND SPECULAR MULTIPATH 


Rodger E. Ziemer 
Electrical Engineering Department 
University of Missouri - RoIIa 
Rolla, Missouri 


ABSTRACT* 

A digital computer simulation of an FC FB loop operating in 
interference and specular multipath is developed. Results from 
several test simulations of the loop operating under various con- 
ditions are presented to demonstrate the usefulness of the model 
and give an indication of the effectiveness of FM spread spectrum 
in combatting multipath and interference. Comparisons are made 
with an ideal discriminator, using normalized mean-square error 
between the output with and without interference present at the in- 
put as the performance measure. These preliminary results indi- 
cate that the FCFB loop is about 100 times better than the discrim- 
inator in suppressing the deleterious effects of RFI and multipath 
under the conditions simulated. The results also are intuitively 
satisfying in that increasing the modulation index of the signal and 
feedback factor of the loop generally lowers the mean-square error 
assuming suitable choices are made which avoid extraneous distor- 
tion effects. 



® age BLAJsIK 
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*This paper was published under X-525-70-329. 
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ERFRAC; A COMPUTER PROGRAM FOR RESEARCH INTO THE 
CORRECTION OF ERRORS IN SATELLITE RADIO TRACKING 
DUE TO ATMOSPHERIC REFRACTIONS 


Bruce H. Morgan 
U. S. Naval Academy 
Annapolis, Maryland 
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* This work was performed in a previous year. 
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ERFRAC: A COMPUTER PROGRAM FOR RESEARCH INTO THE 
CORRECTION OF ERRORS IN SATELLITE RADIO TRACKING 
DUE TO ATMOSPHERIC REFRACTIONS 


Bruce H. Morgan 
U. S. Naval AcacreS^EDING 
Annapolis, Maryland 


blank not mmm 


1. INTRODUCTION 

My previous report on this project appeared in the Final Report of the 1968 ASEE- 
NASA Summer Faculty Fellowship Program (X-520-69-30) on pages 111-132. This re- 
port is cited below as Reference (1). A brief introduction to the problem dealt with is 
given in this section, but the reader is referred to that previous paper for detailed 
derivations and a list of citations to the literature. 

Because of refraction, radio waves transmitted from an earth satellite follow a 
ciu-ved path down through the atmosphere. Consequently, the direction in which they 
are traveling when they reach an observer on the ground is in general different from the 
direction of a straight line from the satellite source to the observer. 

The difference between these directions, is called the angular error at the observer. 
Its magnitude depends in general upon the frequency of the radio signal, the state of the 
atmosphere — troposphere and ionosphere — at the time of observation, and upon the 
angular elevation of the radio source above the observer's horizon. At the Minitrack 
System operating frequency of 136 MHz, is is on the order of one or two minutes of arc 
for an elevation angle of 60°. Since apart from this refraction error, the Minitrack 
operational accuracy is about 20 seconds of arc, it is clearly desirable to apply a re- 
fraction correction to the tracking data. The problem is how to determine this correc- 
tion. Synoptic ionospheric predictions are not adequate for this purpose, particularly in 
view of the frequent occurrence of "TIDs," traveling ionospheric disturbances, which 
cause irregular fluctuation of perhaps one or two minutes of arc in the apparent position 
of a (136 MHz) radio source with a period of half an hour or so. Clearly, then, some 
method of determining the refraction error at the time and place of observation is needed 
if the potential operational accuracy of the Minitrack system is to be attained in practice. 

One possible approach is to observe the directions of arrival of radio waves from 
the same source at two different frequencies. Because the ionosphere is dispersive, these 
two directions will not be the same. The difference between them is a measure of the 
state of the ionosphere along the path from soiu-ce to observer and hopefully enable the 
observer to determine what correction should be applied for refraction. 

My fellowship time this summer (1969) has been devoted to developing a computer 
program (ERFRAC) to trace the path of a radio wave through a model atmosphere which 
includes optionally, a "wedge" of election density at an arbitrary location and orientation, 
Simulating the presence of an ionospheric irregularity. The program is based on the dif- 
ferential equations for the ray path (Euler- Lagrange equations) that derive from Fermat's 
Principle of geometrical optics. This derivation is set forth in Appendix B. By making 
a small change in just two lines of the program, the user can employ either of two 
standard subroutines DRKGS or DHPCG, from the IBM Scientific Subroutine Package, to 
carry out the step by step numerical integration of these differential equations: Sub- 
routine DRKGS uses a Runge-Kutte technique, while DHPCG employs a predictor- 
corrector method. Both are written in double-precision mode, which means that 16 
significant figures are carried in all the calculations. 
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The program as, presently written uses only the simplest ^proximation for the 
relation between electron density and refractivity in the ionsophere. In particular, the 
effect of the earth's magnetic field, which causes double refraction, is not included. A 
better approximation could be incorporated into the existing program and corresponding 
effects, such as Faraday rotation, could be computed. 

The author will be happy to talk with anyone who is interested in makir^ use of 
ERFRAG, in whole or in part. Mail address; Physics Department, U. S. Naval Academy, 
Annapolis, Md. 21402. Telephone: 301-268-7711, extension 393. 


2. GEOMETRY OF THE CALCULATION 

Figiire 1 shows the coordinate system used in the present work. The origin of co- 
ordinates is taken at the center of the earth. The x axis is oriented to pass through the 
point of observation; this direction is referred to in this paper as the observer's vertical. 
The y-axis is taken parallel to the observer's east, and the z-axis, north. This choice 
was made because in the Minitrack system, the direction cosines of the line to a satellite 
source are computed with respect to the observing station's eastward and northward 
base lines. These direction cosines are by convention called L and M, respectively. 

The radio signal originates at the satellite and travels downward through the 
atmosphere to the observer's receiving antenna. For purposes of computation, however, 
it is convenient to consider a hypothetical "reversed ray" which originates at the observer 



Figure 1. Basic Geometry of the Calculation. 


106 



and travels upward. This reversed ray will follow exactly the same path as would a 
downward-traveling ray, since by Fermat's Principle, both will travel that path which 
makes the transmission time an extremum. We simply must start the reversed ray off 
from the observer with an azimuth and elevation angle equal to those that the downward- 
traveling ray has when it reaches him. The path of the reversed ray is described by 
taking x as the independent variable and finding the y and z coordinates of the path as 
functions of x. 


3. SUBROUTINE PERFORMANCE: DRKGS AND DHPCG 

The Euler- Lagrange equations for the ray path are a pair of simultaneovis second- 
order differential equations which can be solved algebraically to give formulas for 
= d^y/dx^ and = d^z/dx^ where the functions y(x) and z(x) are the y and z co- 
ordinates of a point on the ray path as functions of x. 

In order to carry out the integration of these equations to determine y(x) and z(x), 

I decided to make use of one of the subroutines provided in the IBM Scientific Subroutine 
Package. Using the Goddard 360/91 computer, these subroutines are accessible directly 
from the HITS, or CRBE, remote terminals at Goddard. Preliminary calculation had in- 
dicated that double precision (16 significant figure) accuracy might be needed. Two 
suitable subroutines were available in the Package: DRKGS and DHPCG. The first utilizes 
a Runge-Kutta procedure; the second, a predictor- corrector method. 

To check out and compare the operation of these two subroutines, and to become 
familiar with their use, I wrote programs which employed them to integrate a pair of 
second order equations which have simple analytic solutions namely: 




y 


( 12 - 1 ) 


( 12 - 2 ) 


These were integrated from x = 0, to x = 7. For the initial conditions which J 
specified, y(0) = 0.0, y„ (0) = 1, z(0) = 1, and z,(0) = 0.0, the analytic solutions are, of 
course. 


^EXACT = sin X (12-3) 

^EXAcr = X (12-4) 


Therefore, the errors in the computer's calculated values were simply 


and 


Sy = y(x) - sinx. 

(12-5) 

5z = z(x)-cosx 

(12-6) 
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Trials showed that, as might be expected, these errors, at any given value of x 
tended to become smaller as the step size, h = x^+j - X; , used by the integrating sub- 
routine, was made smaller. Also, S y and 8z tended to increase with x, which also would 
be expected. Although each one individually oscillates as a function of x, the sum of their 
absolute values. 


E = (I 8 y| + I 8 z| ), 


exhibited a fairly steady increase with x. 

The results of the trials are summarized in Figure 2. This shows the magnitude 
of the combined error E at x = 6.0 as a function of integration step size h. E is plotted on 
the vertical axis in a quasi- logarithmic manner. On the horizontal axis is plotted the 
subroutine variable IHLF which is the number of times that an initial input parameter is 
cut in half to arrive at the step size used by the integrating subroutine. In all trials shown, 
this parameter was 0.1, so the step-size that was used by the computer is 0.1/2^“'’" . In 
a few trials, the computer gave IHLF one value for some steps and a different value for 
others. The point on the graph for such trials is plotted half-way between the two values 
used. Both DRKGS and DHPCG assign a value to IHLF by comparir^ a "truncation error" 
which they compute at each integration step with an "upper error bound," which is sup- 
plied by the user as an input parameter, called PRMT(4). This parameter affects the 
computation only through its role in setting the value of IHLF. 

The graph. Figure 2, confirms that the error in the result decreases as the integra- 
tion step size is made smaller (IHLF larger) for both DRKGS (circles) and DHPCG 
(triangles). However, both curves flatten out at an error, E, of about 10"’^ . To check 



Figure 2. Combined Error at X = 6.0 vs Step Size Parameter. 
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this behavior, similar curves were plotted in Figure 3, for the error E in the results at 
X = 4.0 radians. (Recall that in figure 2, the error plotted was that at x = 6.0 radians.) 
Again, the two curves each bottom at an E of about 10'’® . The DRKGS curve reaches its 
minimum at a stepsize of about 0.1/2 . 

At both X = 4 and x = 6, then, the result obtained by integration of the differential 
equation is good only to 13 significant figures, although the calculations are all carried 
out with 16 figures. 

Figure 4 shows a similar behavior when DRKGS was used to integrate the equation 
d^y/dx^ = e* with initial conditions y(0) = 1.0 and y (0) = 1.0, for which the analytic 
solution is y = e*. In this case, the error at x = 4.0 reaches a minimum of about 10'“ , 
where the step size used is about 0.1/2®. 

It would be of interest to continue these curves to still smaller step size by making 
additional computer runs. The time available for the project did not permit this. However, 
it appears clear that in all cases, accuracy does not improve when the step-size is re- 
duced beyond a certain point. Whether it actually becomes worse to a significant extent 
is not known, but is quite possible. 

This optimum step size in the cases studies was such that from one thousand to 
four thousand steps filled the range of integration. This suggests that in ERFRAC, where 



Figure 3. 
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Figure 4. Error in y at x = 4.0 vs. Step Size Parameter IHLF, Using Subroutine DRKGS to Integrate 
the equation d^y/dx^ = e from x = 0 to x = 4.0 with initial conditions y(0) = 1.0 and y (0) = 1.0 

the range of interest of the independent variable x is from zero to a few hundred or per- 
haps to one or two thousand kilometers, a step size of one kilometer might be a good 
choice so far as accuracy of the computed ray path is concerned. 


4. COMPARISON OF SUBROUTINES DRKGS AND DHPCG 

A primary objective of the trials described in the previous section was to compare 
the performance of the subroutines DRKGS and DHPCG in order to decide which one to 
use in ERFRAC. Possible criteria that could be used in making this comparison include 
the degree of accuracy which is attainable; the amount of computer time required to 
produce a given degree of accuracy in the output; and the dependability or predictability 
of achieving a given degree of accxmacy. The trials indicate that DHPCG and DRKGS are 
about equal under the first criterion, DHPCG is definitely superior under the second , 
criterion, and DRKGS is perhaps somewhat superior under the third. My conclusion 
from these trials was to utilize DHPCG. 
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With regard to the first criterion, Figures 2 and 3 show that, while the error in the 
output of DHPCG using a given integration step-size was less than that of DRKGS, both 
DRKGS and DHPCG appear to have about the same ultimate limit of accuracy. This limit 
for each was at a "combined error," E, equal to about 10‘*^ and occurred for both at about 
the same step-size, namely h = 0.1/2 . 

With regard to the second criterion, we note first from Figure 2 that for a given 
step-size, the error in the output of DRKGS is from five to ten times as large as that of 
DHPCG, at least down to a step- size where both subroutines approach their ultimate limit 
of accuracy. (This limit, as pointed out above, is about the same for each.) Secondly, 
calculations in Appendix G, using the computer time charged to the trial run jobs, indicate 
that, again for a given step size, the time required by DHPCG to complete a run over a 
given range of x is about one third less than that required by DRKGS. So for a given step 
size, DHPCG is both more accurate and takes less time. 

This comparison can be made somewhat more quantitative; by the following 
reasoning: Since the time required by DHPCG at a given step-size is 1/3 less than that 
required by DRKGS, DHPCG could use a step-size 1/3 smaller than DRKGS does and 
finish in the same time. But from Figure 2, a step-size 1/2 as large leads to an error 
of only about 1/12 as much. A reduction in step-size of 1/3 should therefore reduce 
error by a factor of about 2/3 of 12, or 8. Since the error of DHPCG at the original 
step-size was already only 1/5 to 1/10 of that of DRKGS, the conclusion is that in a 
given time, DHPCG could produce results with errors only some 1/40 to 1/80 as great 
as DRKGS. 

Finally, with regard to the third criterion, it appears from Figure 2 that the trend 
toward smaller error with smaller step size is somewhat more regular for DRKGS than 
for DHPCG. In other words, the data points plotted there for DRKGS fall more nearly 
on a straight line than do those for DHPCG, until the ultimate limit of accuracy is ap- 
proached. This indicates that the accuracy of DRKGS is somewhat more predictable 
than that of DHPCG. I believe, however, that the difference in this respect shown by 
Figure 2 is not large and is not very important, and that therefore this factor is not 
entitled to much weight. 


5. THE BASIC ALGORITHM OF ERFRAC 

As mentioned above, the basis for ray tracking used in this project is a pair of 
second order differential equations for the ray path derived in Appendix B of Reference 
(1). These equations relate the curvature of the path, expressed in rectangular coordinates, 
to the slope of the path and the magnitude and space derivatives of the local radio refrac- 
tive indeix, y-. The subroutines DRKGS and DHPCG are each designed to integrate a sys- 
tem of first-order differential equations. But two second-order equations can easily be 
written as a system of four first-order ones by a well-known reformulation. 

First, since the refractive index and its derivatives are functions of position; 
we can write our two second order equations as follows: 


y/d^ X = g (x. y, z. y^, 


(18-1) 


d^ z/d^ X = h (x, y, z, y^, z^) 


(18-2) 
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In order to conform with the variable names used in DRKGS and DHPCG, we now 
define: 


y = Vj and z = y^, 


(18-3) 


and we also introduce: 


y, = d y/d X = y^ and = d z/d x = y_^ 


(18-4) 


Then we can write 


d y^/ d X = d y/d x = y^ 


(18-5) 


d y/d X = d2 y/d x2 = g (x, y^, y^, y^, y^) 


(18-6) 


d y /d X = d z/d x = y^ 


(18-7) 


d y/d x = d^ z/d x^ = h (x, y^, y^. y^, y^) 


(18-8) 


These last four equations are a set of first order equations with independent vari- 
able X and dependent variables Yj , , Yj , and y^ . They constitute the basis of the actual 

computer program. 


6. PROGRAM LOGIC IN ERFRAC 

The purpose of this section is to provide an explanation of the logic in ERFRAC. A 
flow chart of the program is given in Fig. 5 and a listing of the source deck in Fig. 6. The 
reader may also refer to the list of variable names and other program terminology given 
in Appendix A. 

The explanation in this section is keyed to the 'TSN' numbers which appear in the 
left margin of the source deck listing in Fig. 6. 

(Note: The argument lists for both of the user-provided subroutines; FCT and 
OUTP, are specified by the IBM subroutine DHPCG. My program ERFRAC, as written, 
uses variables in these lists to pass initialization data from the main program to the sub- 
routines and from one subroutine to another. This could alternatively be done by the use 
of COMMON storage.) 

ISN 

Main Program Comment 

5a, 5b Comment mark avoids execution of READ statement when data is speci- 

fied in the program, as it is for this run (in lines 9 to 13). 
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ISN 

Main Program 
6-8 

14, 15 

16 

17 

18, 19 
20 
21 

22-24 

25 

26 
27 


FCT 

7 

8-20 


Comment 

Set up DO loops for this particular run, to sample effect of various 
values of HWGB (line 23), AZWG (line 24), and FREQ (line 22) 

These steps were provided so that the start and end of the range of x 
over which the radio ray path is to be traced can be specified on input 
in terms of altitude above the earth's surface rather than in terms of 
distance from the earth's center. 

The value of RO is transferred to the subroutines FCT and OUTP as X. 

The value of PRMT(l) =XO is transferred to subroutine FCT as Y(2). 

The reversed ray path starts at the observing station which by choice 
of the coordinate system lies on the x-axis, so y = Y(l) = O and z = 

Y(3) = O. 

PRMT(6) = 0 causes subroutine OUTP to perform a sequence of 
initialization steps, one of which is to change PRMT(6) to unity so 
that they are shipped over on subsequent calls. 

DERY(l) serves the same purpose for the subroutine FCT that PRMT(6) 
does for OUTP. (See comment just above.) A large negative value is 
used for initialization because later on, DERY(l), is the value of dy/dx 
along the ray path and this is unlikely to become large and negative. 

See 6-8 above. 

At this point the program takes up initialization procedure in subroutine 
FCT, lines 3 to 66, as described below. 

Now we carry out initialization in subroutine OUTP, lines 2 to 38, as 
described below. 

Subroutine DHPCG is called to carry out the first step of tracing the 
path of the reversed ray. It calls on subroutine FCT to calculate the 
local values of the space derivatives of the path, and hands the co- 
ordinates of the new point on the path to subroutine OUTP which 
processes them and then returns control to DHPCG. This cycle, 
DHPCG to FCT to DHPCG to OUTP to DHPCG, is repeated until the 
end value of x, PRMT(2), as specified on input, is reached. Then 
control returns to the MAIN program where, in the version listed in 
Fig. 6, a new set of input parameters is fixed and the entire program 
is repeated, until the DO loops of the MAIN program are completed. 


Checks whether this is the initial call. If not, the initialization procedure 
in lines 7A to 66 is skipped. 

On the run shown, data was specified in the program rather than being 
read in from data cards. 
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Figure 5(a) 
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FCT 



END 

Figure 5(b) and (c) 
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(FCT) 


23-29 

b 


30 

31 


Comment 


The value of RO ane XO are transferred from MAIN to FCT as X and 
Y(z), respectively. 

If there is no wedge of electron density (MEWGMX = O), then HWGB 
is set to a large positive value. This causes the program at line 80 
below to skip the electron wedge calculations and go directly to line 
111. TWGO, SLOPWG, and AZWG are also given large values so that 
asterisks will be printed for the values of these quantities in the out- 
put data heading. 

Write main heading for output data. 

Set numerical multiples for converting angles from degrees to radians. 


32-35 

36 

37-39 

40-41 

42 

44-50 


Convert ALTOBS and AZOBS from degrees to radians and store in 
DERY(2) and DERY(3) for transfer to OUTP. 

S e (ALTOBS) in DERY(4) for transfer to OUTP. 

SINAZO, and COSAZO. 

Com. - SINAXO and COTALT - COSAXO and store in 

Y(2) ana n4; for transfer to OUTP. 

If there is no electron wedge the initialization of FCT is complete and 
the program returns to MAIN. 

Compute the radius vector of the point where the projection of the 
reversed ray intersects the base of the electron wedge. (See Appendix D). 
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Comment 


ISN 

(FCT) 


51-55 


56-58 

59-62 


63-65 


66 


67-69 


70 


71 

72-74 

75 

76-73 


80 

82-87 


87A 


88-90 


93 


If the value held by the computer for RWGB is less than XWGB, then to 
avoid an error signal at 58, the values of COSTHT and SINTHT are set 
rather than computed. (NOTE: LT. in 51 should be LE.) 

Computation of COSTHT and SINTHT . 

Convert . AZWG from degrees to radians and compute SINAZW and 
COSAZW. 

Compute AWG, BWG, CWG. (See Appendix E.) 

Return to MAIN; initialization of FCT is completed. 

Use output from DHPCG to compute altitude, H, of current point on ray 
path. 

Compute tropospheric radio refractivity. (Reference: Bean, B.R., and 
Dutton, E. J., Radio Meteorology (NBS Monograph 92) GPO 1966, 
pp. 65-67). 

Conipute gradient of tropospheric refractive index. 

Compute electron density, NE, of Chapman layer. (Reference: Davies, 
K., Ionospheric Radio Propagation (NBS Monograph 80) GPO, 1965, 
page 17.) 

Compute refractive index due to Chapman layer. 

Compute gradient of the refractivity due to the Chapman layer and the 
troposphere, superimposed. 

If the ray path has not yet reached the altitude of the base of the elec- 
tron we(%e, then skip to 111. 

Compute the height of the ciirrent point relative to the base plane of 
the wedge (HWG) and the local thickness of the wedge (TWG) 

Call for a print-out, used in debugging the program. 

If HWG is equal to or greater than TWG, then we are through the wedge. 
In that case, we set HWGB at a large value so on subsequent calls the IF 
statement at 80 will cause all computations involving the wedge to be 
skipped over, and proceed to 111. 

If HWG is less than zero, the reversed ray has not yet entered the elec- 
tron wedge, and we skip to 111. (It is possible, though not likely, for 
HWG to be less than zero even though H is greater than HWGB because 
the base of the wedge is taken as plane so that its altitude above the 
earth’s surface is not constant. HWGB is its minimum altitutde. See 
Fig. D-1). 
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ISN 

(FCT) 

95-104 

104A 

105-108 

109-110 

111-112 

113-114 

114A-115 

116-135 

OUTP 

6-10 

11,13 

14-22 

23-24 

24A-29 

30, 31 
32, 33 


Comment 


Compute the electron density, NEWG, and the components of the gradient 
of the refractive index due to the electron wedge at the current point. 

Statements used in debugging the program, to produce a line of print-out 
at this point each time subroutine FCT is called by DHPCG. 

The total electron density and the gradient of the refractive index at the 
ciirrent point are the respective sums of the values due to the Chapman 
layer (computed at 74 and 77-79) and those due to the electron wedge. 

The local refractivity is the sum of that due to the total electron density 
and that associated with the troposphere (line 70). 

If electron wedge calculations have been skipped over, then the local 
refractivity is just the sum of that due to the Chapman layer and that 
due to the troposphere (line 70). 

In any case, the refractive index follows from the refractivity, by defini- 
tion of the latter. 

Statements used for debugging to produce a line of print-out at this point 
each time subroutine FCT was called by DHPCG. 

Computation of DERY for return to the calling subroutine, DHPCG. The 
formulas used are derived in Appendix B of Reference (1) cited in the 
Introduction (Goddard Space Flight Center, X-520-69-30). For convenience, 
that Appendix is reproduced as Appendix B of the present paper. Note 
However that the independent variable used in ERFRAC is x, not z, so 
that in the formulas of Appendix B, z must be changed to x, x to y, and 
y to z. 


Format statements; will be mentioned where the read and write state- 
ments appear below. 

If this is not the initial call of OUTP, the initialization proceudre, lines 
14-38, is skipped over. 

Values of variables are set. 

LOBS and MOBS are computed; see App. F. 

On the run listed, the values for DERY and for DXPRNT are specified 
in the source deck rather than read in from data cards. The appropri- 
ateness of the values used here for DERY, the error weights used by 
DHPCG, is uncertain and should be checked. The value of DXPRNT 
was chosen so as to obtain at least one output data point in the electron 
wedge, whose thickness was set at 10 km. 

Initial values of H and XMXO are computed for first line of print-out. 
Print column headings and initial line of data. 
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INS 

(OUTP) 


Comment 


34-36 Set the three print control indexes to unity. 

37 Set PRNTFX at a positive value which is less than the smallest step 
size that may be used in DHPCG, namely PRMT(3)/2'“ . 

38 Return from OUTP to MAIN program following initialization of OUTP. 

39 On each call subsequent to the first, execution begins here (Fortran 
label 70). 

40 Compute chaise in X, since the last call, DX. 

41 Return to DHPCG if DX is zero. If DHPCG, when it is first call, merely 
hands the initial values of the coordinates of OUTP, then OUTP returns 
to DHPCG for the first step of integration. 

43-45 Computes the increment of path length, DS, for the current step of 

integration. (This is not utilized in the present run, but could be ac- 
cumulated and printed-out if desired. Also, if in subroutine FCT, 
the current value of the refractive index MU were stored in COMMON, 
this could be accessed by OUTP to multiply DS, giving the increment 
in optical path length. 

46-48 Store current values of coordinates of point on ray path. 

49 Compute XMXO 

50 Retxu-n to DHPCG if it is not time for a line of print-out. 

52 Increase line print control index by one. 

53-57 Print a blank line after every fifth line of output. 

58-63 Provides pagination, with column headings on each page. Note: In 

order to obtain an integral number of groups of five lines each on a 
page, "IPRNT" in line 58 should be replaced by ’TPRNT-2." 

64-65 Compute R and H for current point. 

66-69 Compute displacement of current point from point at same value of 

X on a straight line projection of the ray as observed at the observing 
station. 

70-71 Compute RANGE, the distance from observing station to current point. 

72-77 If the value held in storage for XMXO is equal to or greater than that 

for RANGE, then the true elevation angle of the current point is directly 
set to 90°, thus avoiding the generation of an error signal in an attempt 
to calculate it, at 85. Also AZS is set to a value which will cause aster- 
isks to be printed for azimuth error in the print-out data. 

78-79 Compute the true elevation angle of the current point, and its distance 

from the X-axis, RPOLS. 
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INS 

(OUTP) 

80-84 


85 

87-91 


92-95 


96-108 


109-110 


Comment 


If the computer's value of y is less than its value of RPOLS, then we 
avoid an error signal by setting the azimuth of the current point to 90° 
without trying to compute it. 

Compute the azimuth angle of the current point. 

Compute the difference between the observed and actual elevation angle 
of the current point in radians, DEVALT, milliradians, DALTMR, and 
seconds of arc, DALTS; and the difference between the observed and 
actual azimuth in radians, DEVAZ, and seconds of arc, DAZS. 

Compute the true direction cosines for the current point and their dif- 
ferences from those of the ray direction as observed at the observing 
station. 

Compute the angle between the tangent to the ray path at the current point 
and the direction of the ray as observed at the observing station, in 
seconds of arc. 

Write a line of output data and retimn to DHPCG. 


7, RESULTS 

Successful complete runs of ERFRAC, tracing a radio-ray up and through electron 
wedges, were achieved on November 10, 1969. Figures 7 through 11 show results obtained 
on a run for which the parameters were as follows. (See list of symbols in Appendix A 
for definitions.) 

Radio frequency, FREQ = 136 MHz, 

Observed Elevation ai^le of ray, ALTOBS = 50 degrees 

Observed azimuth of ray, AZOBS = North 80 degrees East 

Chapman Layer Parameters; 

Maximum electron density, NEMAX = 2 x 10*^ n’^ 

Height of maximum density, HEMAX = 300 km. 

Scale Height of Chapman layer, HCMP = 100 km 
Secant of sun's zenith angle, SECSUN = 1,0 

Electron Wedge Parameters; 

Height of base of wedge, HWGB = 250 km 
Thickness of wedge where ray enters, TWGO = 10 km 
Tangent of vertex angle of wedge, SLOPWG = 0.1 
Azimuth of wedge gradient, AZWG = 0 degrees 
Maximum electron density, NEWGMX = 1 x 10^^ m"^ . 

Figure 7 shows to an exaggerated scale the way in which the radio ray is caused to 
curve by the varying refractive index of the atmosphere. As indicated, the graph shows 
only the eastward component of the horizontal displacement of the ray from a straight line 
projection, but this is nearly the total horizontal displacement since the observed azimuth 
angle of the ray was set to 80°, so the northward component is small (less than 1/5 of the 
eastward component, according to the print-out). 
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Figure 7. Radio Ray Path, Showing Eastward Component of Horizontal Displacement from 
Straight Like Projection of the Observed Ray Direction, at the Greatly Exaggerated Scale 
Indicated, 

Figure 8 is another plot of this eastward horizontal displacement, and Fig. 9 shows 
a magnified portion of Fig. 8, revealing more clearly the effect on the ray path of the 
assiuned electron wedge at 250 km altitude. 

Figure 10 shows the elevation angle error (the difference between the apparent and 
actual elevation of the radio source) as a function of the altitude of the source. The effect 
of tropospheric refraction is clearly apparent, producing an angle error of some 40 
seconds of arc for a source only 20 km up. On the other hand, the effect of the electron 
wedge at 250 km altitude is hardly discernible. Its presence shows up in Figure 11, how- 
ever, where the total or cumulative bending is plotted vs. altitude for a ray which starts 
downward from an altitude (XMXO) of 400 km. It gradually curves upward as it travels 
down into the Chapman layer until it reaches the altitude of maximum electron density 
(300 km) and then it begins curving downward again. The electron wec^e causes a sharp 
bend upward and then a sharp downward bend which is not quite as great so that the 
cumulative bending below the wedge is more positive than it would have been without the 
wedge. Consequently, for this wedge orientation (AZWG = 0 degrees), the magnitude of 
the elevation angle error at the observing station is less than it would have been without 
the wec^e. An orientation with AZWG = 180 degrees presvimably would cause maximum 
ai^le error. (The angle error with AZWG = zero was 453 seconds of arc; with AZWG = 90°, 
it was 470 seconds, and with no wedge at all 464 seconds.) 
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HORIZONTAL DISPLACEMENT OF RAY FROM 
STRAIGHT LINE PROJECTION (meters) 


Figure 8. Eastward Displacement of Ray Path from a Straight 
Line Projection of the Observed Ray Direction 


The cumulative bending becomes negative at altitudes below 200 km, indicating that 
the downward slant of the ray path is steeper than it was when it started at 400 km. (The 
slant ai^le which correspond with this plot is the angle of the ray relative to the observ- 
ing stations, horizontal rather than local horizontal.) 

The integrating subroutine, DHPCG,used the step-size specified on input, PRMT(3) 

= 1 km throughout the run except at the start and in the electron wedge. At the first 
print-out point, 2 km altitude, it was using step-size 1/2^ = .125 km, and increased this 
to the nominal 1 km at 44 km altitude. When it encountered the large electron density 
gradients in the electron wedge, it decreased the step- size, from 1 km at altitude 249 km, 
to 1/2^ = .031 km at 251 km, returning to 1 km step-size at 266 km, which is above the 
wedge. (DHPCG reduces the step-size if an internally- computed measure of "truncation 
error" at a single step exceeds an "upper bound" which is specified as PRMT(4) on 
input.) 
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STATION 

ALTITUDE 

(km) 



HORIZONTAL DISPLACEMENT OF RAY FROM 
STRAIGHT LINE PROJECTION (meters) 


Figure 9, Magnified Portion of Figure 8, Showing 
Effect of Electron Wedge 

Results for the elevation angle error for a variety of circumstances are shown in 
the table below. The last three columns show that when the angle error ascribable to 
tropospheric refraction (about 50 arc seconds), is subtracted, the remaining angle error 
at operatii^ frequency 68 MHz is very nearly just four times that at 136 MHz under all 
conditions, and is constant at least to slide- rule accuracy, for a given satellite altitude. 
From this, it appears that if the apparent elevation were observed at both frequencies, 
the elevation angle error due to the ionosphere could be computed to within one percent 
which is about 4 arc seconds in the present case, by taking the difference between the 
two observed elevation angles and dividing by 3.02 or thereabouts, the exact value de- 
pending on satellite altitude. The result would then be added to the elevation observed 
at the lower of the two observed frequencies to get the true elevation. 

Of course, careful theoretical analysis is needed before the adoption of such a 
procedure could be seriously proposed. Such analysis is not attempted in this paper. 

For reference, the errors in the direction cosines L and M (these are the cosines 
of the angles measured from station east and from station north, respectively), which 
were also computed and printed out, are tabulated below. The values tabulated are 10^ 
times the actual error. For example, in the first line, the error in L at 68 MHz is 612 x 
lO'® = .00612, and the error in M at 136 MHz is 29,2 x lO"® = .000292. 
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IONOSPHERIC REFRACTION- 


Figure 11. Cumulative Bending of Ray vs. Altitude, 
for Source at 400 km Altitude. 





Satellite 

Altitude 

(km) 



Azimuth 
Angle of 
Electron 
Wedge 


Elevation Angle 
Error, e (arc sec) 


“ ^TROP 

(arc seconds) 



68 MHz 

136 MHz 

68 MHz 

136 MHz 


453 

— 

403 


458 


408 


470 

1699 

420 


443 

1583 

393 

1667 

451 

1617 

401 

1744 

470 

1694 

420 

1383 

382 

1333 

332 

1424 

392 

1374 

342 

1514 

414 

1464 

364 

1166 

328 

1116 

278' 

1211 

339 

1161 

289 

1310 

364 

1260 

314 


Ratio of 
last 2 cols. 



Satellite 



(Error in Direction 
Cosines) x 10® 


Altitude 

(XMXO) 

AZOBS 

Error in L 

Error in M 



68 MHz 

136 MHz 

68 MHz 

1 

136 MHz 

400 

0? 

612 

166 

108 

29.2 


45° 

617 

167 

128 

34.0 


90° 

633 

171 

137 

36.2 

600 

0° 

595 

162 

105 

28.5 


45° 

603 

164 

136 

36.2 


90° 

629 

170 

151 

39.8 

800 

0° 

504 

140 

88.9 

24.6 


45° 

513 

142 

126 

33.7 


90° 

544 

149 

143 

37.9 

1000 

0° 

425 

120 

75.0 

21.2 


45° 

435 

122 

115 

31.1 


90° 

469 

131 

134 

35.8 
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APPENDIX A 
LIST OF VARIABLES 


Note; All distances are measured in meters. 


A through H 

In termediate variables used in computing the values of DERY in 
the subroutine FCT when called upon by DHPCG or DRKGS. 

A, Al, B 

Intermediate variables used in computing the space derivatives 
of the refractive index in the electron wedge; in FCT. 

ALTOBS, AZOBS 

The elevation angle and the azimuth angle, measured from north 
toward east, of the radio ray as observed at the observing station; 
in OUTP. ALTOBS and AZOBS also appears in FCT, first in 
degrees as input parameters, and then in radians. Their values 
are passed to OUTP as DERY(2) and DERY(3), respectively. See 
Figure C-1. 

ALTS, AZS 

The elevation angle and the azimuth angle, measured from north 
toward east, of the radiiis vector from the observing station to the 
current point on the ray path; in OUTP. 

AUX 

A two-dimensional storage array used internally by subroutines 
DRKGS and DHPCG. For ERFRAC, its dimensions are (8,4) if 
DRKGS is used, and (16,4) if DHPCG is used. 

AWG, BWG, CWQ 

Direction cosines of the electron wedge gradient, in FCT. 

AZOBS 

See ALTOBS. 

AZS 

See ALTS. 

AZWG 

An input parameter (in degrees) which specifies the azimuth angle 
of the density gradient of the electron wedge; in FCT. (Its value 
is initially transferred from the MAIN program as DERY(4).) See 
App. E. 

B 

See A. 

BNDG 

The angle, in radians and then in seconds of arc, between the 
tangent to the ray path and the direction of the ray as observed at 
the observing station; in OUTP. 

BWG 

See AWG. 

C 

See A 

CHPTRM 

An intermediate parameter which appears in the formula for the 
electron density in a Chapman layer; in FCT; CHPTRM = SECSUN 
* DEXP (-ZCHP). 

COSAZO 

The cosine of AZOBS; in FCT. 

COSAZW 

The cosine of AXWG; in FCT. 
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COSBDG 


Cosine of BNDG 


COSTHT Cosine of the angle between the X axis and the line from the earth's 

center to the electron wedge base point P; in FCT. See Fig. D-1. 

COTALT Cotangent of the elevation angle of the radio ray as received at the 

observing station; in FCT. 

CWG See AWG. 

D See A. 

DALTMR This is DEVALT, in milliradians. 

DALTS This is DEVALT in seconds of arc; in OUTP. 

DAZS This is DEVAZ in seconds of arc; in OUTP. 

DENOM An intermediate variable in the evaluation of DERY(2) and DERY(4); 

in FCT. 

DERY A four-element one dimensional array. At the time of the first call 

of the integrating subroutine (DRKGS or DHPCG), it is the set of 
"error weights" used by the subroutine in computing the truncation 
error at each integration step. Subsequently, its elements are the 
derivatives of the independent variables, Y, as computed by sub- 
routine FCT. 

Before the integrating subroutine is ever called, DERY is used in 
initialization procedures, as follows: 

DERY(l) is given the value -10’® as a way for FCT to determine 
whether it is being called for the first time. 

DERY(2) is first used to transfer the value of FREQ from the 
main program to FCT, when the main program contains a DO Loop 
to repeat runs with different values for FREQ. Then it is used to 
transfer the value of ALTOBS from FCT to OUTP. 

DERY(3) is first used to transfer the value of HWGB from the 
main program to FCT, when the main program contains a DO Loop 
to repeat runs with different values for HWGB. Then it is used to 
transfer the value of AZOBS from FCT to OUTP. 

DERY(4) is first used to transfer the value of AZWG from the main 
program to FCT, when the main program contains a DO Loop at 
repeat runs with different values for AZWG. Then it is used to 
transfer the value of the sine of ALTOBS from FCT to OUTP. This 
transfer is not absolutely necessary, but the sine of ALTOBS is 
needed for initialization computations in both FCT and OUTP and 
the transfer merely avoids having to calculate it twice. 

DEVALT Difference between the elevation angle of the radio ray as observed 

at the observing station and the true elevation angle of the current 
point on the actual ray path, in radians; in OUTP. A positive value 
means the observed elevation is greater than the true elevation. 
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DEVAZ 


DEVL, DEVM 

DEVM 

DGTORD 

DHPCG 

DRKGS 

DS 

DX, DY, DZ 
DXPRNT 

DXWG, DYWG, 
DZWG 

DY 

DYDXO, DZDXO 

DYWG 

DY2NUM 

DY4NUM 

DZ 

DZDXO 

DZWG 


Difference between the azimuth angle of the radio ray as observed 
at the observing station and the azimuth of the current point on the 
actual ray path, both measimed from north around toward east; in 
OUTP. A positive value means the observed azimuth is greater 
than the true azimuth. 

Errors in the direction cosines of the radio source as measured 
from the y (east) axis and from the z (north) axis, respectively; in 
OUTP. A positive value means the true direction cosine is greater 
than the one observed. 

See DEVL. 

Numerical Constant which, as a multiplier, converts arc degrees 
into radians; in FCT. 

A double-precision Fortran- language subroutine listed and explained 
in the IBM Scientific Subroutine Package (1 968) , beginning on page 
337. It uses Hamming's modified predictor-corrector method to 
numerically integrate a system of first-order ordinary differential 
equations with given initial values for the dependent variables. 

A double-precision Fortran- language subroutine listed and explained 
in the IBM Scientific Subroutine Package (1 968) beginning on page 
333. It uses the Runge-Kutta method to numerically integrate a 
system of first-order ordinary differential equations with given 
initial values for the dependent variables. 

Increment of length aloi^ the ray path, in OUTP, 

Rectangular components of DS; in OUTP. 

An input parameter which specifies the increment of the independent 
variable, x, at which output is to be printed; in OUTP. 

Rectangular Coordinates of point in electron wedge relative to the 
wedge base point; rectangular components of the vector R in 
Fig. H-1; used in FCT. ® 

See DX. 

Rate of change of y and of z, respectively, with respect to x, along 
the direction of the reversed ray path as observed at the observing 
station; in OUTP. The initial values of Y(2) and Y(4), respectively. 

See DXWG. 

An intermediate variable in the evaluation of DERY(2); in FCT. 

An intermediate variable in the computation of DERY(4); in FCT. 

See DX. 

See DYDXO. 

See DXWG. 
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E 

ERFRAC 

F 

FCT 

FREQ 

G 

H 

H 

HCHP 

HEMAX 

HWG 

HWGB 

IHLF 

IP AGE 
IPRNT 
ISP ACE 

JAZWG 

JFREQ 


See A. 

The name used by the author for the entire computer program 
which is the subject of this document. See flow chart, Fig. 5. 

See A. 

The name of the subroutine called in both DRKGS and DHPCG to 
calculate the first derivatives with respect to X of the four de- 
pendent variable in the array Y. 

The frequency of the radio signal being transmitted, in megahertz; 
in FCT. Its value is transferred from the MAIN program to FCT 
as DERY(2). 

See A. 

See A. 

Height measured vertically (radially from center of earth) above 
the surface of the earth, H = R - RO; in OUTP. 

Scale height in the formula for the electron density in a Chapman 
layer; in FCT. 

Parameter in the Chapman layer formula for electron density 
which represents the height above the earth's sxirface at which the 
electron density is maximum ; in FCT. 

Height above the base plane of the electron wedge; in FCT. See 
Fig. H-1. 

Height above the base of the election wedge perpendicular to this 
base; in FCT. See Fig. H-1. 

This integer gives the number of times that the initially specified 
integration step size - PRMT(3) - has been cut in half by DRKGS 
or DHPCG before carrying out the step just completed. It is 
handed to subroutine OUTP, which can print it if desired. 

A print control index used to control the paging; in OUTP. 

Line printing control index; in OUTP. 

A print control index used to provide a blank line after every five 
lines of print, in OUTP. 

The integer variable index of a DO Loop used to repeat computer 
runs with different values for AZWG in one mode of executing 
ERFRAC. 

The integer variable index of a DO Loop used to repeat computer 
runs with different values for FREQ, in one mode of executing 
ERFRAC. 
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JHWGB 

L 

LOBS, MOBS, 
NOBS 

LPATH, MPATH, 
NPATH 

LSQ 

LTRU, MTRU 

MOBS 

MPATH 

MTRU 

MU 

MUOVG 

MUR 

MUX, MUY, MUZ 

MUXWG, MUYWG, 
MUZWG 

MUY 

MUYWG 

MUZ 

MUZWG 

NDIM 

NE 


The integer variable index for a DO Loop used to repeat computer 
runs with different values for HWGB, in one mode of executing 
ERFRAC. 

An intermediate variable used in computii^ XWGB; in FCT. 

Direction cosines of the reversed radio ray as observed at the ob- 
serving station; measured from the y(east, or y(l)) axis, the z 
(north, or Y(3)) axis, and the x (vertical) axis, respectively; in 
OUTP. 

Direction cosines of the currently computed increment of ray path, 
DS, in OUTP, as measured from the y (east), z (north), and x (station 
vertical) axes, respectively. 

L * L, 

The direction cosines of the radius vector from the observing 
station to the current point on the ray path relative to the y (east) 
and z (north) axes, respectively. 

See LOBS. 

See LPATH. 

See LTRU. 

The atmospheric radio refractive index, in FCT. 

MU/G; in FCT. 

Gradient of the atmospheric radio refractive index. It is computed 
as the sum of the gradient of the tropospheric refractive index, 
which is evaluated first, and the ionospheric refractive index, 
from the Chapman layer formula. 

The partial space derivatives of the refractive index associated 
with the troposphere and the Chapman layer; in FCT. 

Spatial derivatives of the refractive index associated with the 
electron index associated with the electron wedge; in FCT. 

See MUX. 

See MUXWG. 

See MUX. 

See MUXWG. 

The number of first-order differential equations in the system to 
be integrated; it is 4 in ERFRAC, and is so specified in the MAIN 
program. 

Electron density due to the Chapman layer plus density due to elec- 
tron wedge, if any; in FCT. 
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NEMAX 

NEWG 

NEWGMX 

NMU 


NMUCHP 

NOBS 

NPATH 

OUTP 

PRMT 


PRMT(l) 

PRMT (2) 

PRMT (3) 
PRMT (4) 


The maximum value of the electron density associated with the 
Chapman layer - an input parameter in units of 10*® electrons/ 
m^; in FCT. 

Electron density in the electron wec%e, in units of lO^o electrons/ 
m^; in FCT. 

An input parameter which specifies the maximum electron density 
associated with the electron wedge in units of 10*® electrons/m®; 
in FCT. See Appendix H, Equation H3-1. 

Radio refr activity of the troposphere; in FCT. It is related to the 
refractive index, MU, by the following formula: 

NMU = (M U - 1) * 10® 


See App. I. 

Radio refr activity due to the Chapman layer; in FCT. See App. I. 

See LOBS. 

See L PATH. 

The name of the external subroutine called in both DRKGS and 
DHPCG to receive and utilize the computed results at each inte- 
gration step. 

A six-element one dimensional array, as listed below. Only the 
first five elements of PRMT are needed or used by the subroutines 
DRKGS and DHPCG. The sixth element is used in ERFRAC as an 
initialization index for subroutine OUTP. Additional elements can 
be added if desired by simply dimensioning PRMT greater than six. 
Such elements could be used to transfer additional numbers from 
the main program to subroutine OUTP . 

This is the value of x at which integration of the reversal ray path 
begins - the x-coordinate of the station which observes the ray. It 
is also called XO. To facilitate input, the height of the station rela- 
tive to the earth, of radius RO, was chosen as the quantity to read 
in. RO is then added by the computer to obtain PRMT(l). 

This is the value of x at which integration of the reversed ray path 
terminates. It is called XEND in the integrating subroutines. To 
facilitate input, the quantity XEND minus RO was chosen as the 
quantity to read in. RO is then added by the computer to obtain 
PRMT (2). 

The initial value of the integration step size, set by the user. 

The "upper error bound," used in each of the integrating subroutines 
DRKGS and DHPCG in setting the value of the step size parameter 
IHLF. See discussion in Sec. 3 of this paper. 
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PRMT(5) 

PRMT(6) 

PRNTFX 

QWG 

R 

RO 

RANGE 

RANGSQ 

RPOLS 

RWGB, XWGB, 
YWGB, ZWGB 

RWGBSQ 
SEC SUN 

SINAZO 

SINAZW 

SINTHT 

SLOPWG 


DRKGS and DHPCG initialize this parameter to zero. The user can 
stop the integration at any output point by changing it to non-zero in 
subroutine OUTP. 

PRMT(6) is initialized by the main program to the value O. As a 
way of informing subroutine OUTP whether or not OUTP is being 
called for the first time, for initialization purposes. At the first 
call, of OUTP, PRMT(6) is set to unity. 

This parameter was introduced in a line print control instruction to 
get output printed at the discrete values of the independent variable 
X which are closed to the specified values, IPRNT * DXPRNT, 
rather than at the x used in the next succeeding step of integration; 
in OUTP. 

Dimensionaless Parameter which varies from -1, at base of elec- 
tron wedge, to +1 at its top, along a line perpendicular to the base 
plane of the wedge; in FCT. See App. H. 

Distance from the center of the earth to the point at which the cur- 
rent integration step is being performed; in FCT and in OUTP. 

The radius of the earth, which is assumed spherical in ERFRAC. 

RO is the datum for calculating the altitude in the formulas for 
tropospheric and Chapman- layer (ionospheric) refr activity; in 
FCT. The variable X is used to transfer the value of RO from the 
MAIN program to FCT and OUTP, 

Straight-line distance from observing station to cimrent point on 
ray path. 

RANGE - RANGE. 

The component of the RANGE perpendicular to the x (station 
vertical) axis; in OUTP. ("POLS" signifies "polar coordinate, rela- 
tive to observing station.") 

Magnitude and rectar^lar components of the radius vector from the 
center of the earth to the electron wedge base point; in FCT. See 
Fig. H-1, 

RWGB * RWGB. 

Secant of the zenith angle of the sun; appears in the formula for 
the electron density in a Chapman layer, in FCT. 

Sine of AZOBS. 

Sine of AZW. 

Sine of the angle between the x axis and the line from the earth 
center to the electron wedge base point, P, in FCT. See Fig. D-1. 

Tai^ent of the vertex angle of the electron wedge; in FCT. 
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SWG 


Distance parameter within the electron wedge, in FCT. See Fig. 
H-1 and App. H. 

TWG Thickness of the electron wedge, measured perpendicular to the 

base plane of the wec^e, at the point of calculation; in FCT, See 
Fig. H-1. 

TWGO Vertical thickness of the electron wedge at the wec^e base point; 

in FCT. See Fig. H-1. 


WGPl, WGP2 
X, X 


XO or Xg 

XMRO 

XMXO 

XOLD, YOLD, 
ZOLD 


Intermediate parameters used in computing the value and space 
derivatives of the refractive index associated with the electron 
wecfee, in FCT. 

The independent variable in the ray tracing computation: the x- 
coordinate of position, measured relative to a rectangular refer- 
ence frame having its origin at the center of the earth, with its 
x-axis passing through the station observing the ray, X is the 
variable name used for x in ERFRAC. Also used initially to trans- 
fer the value of RO from the MAIN program to FCT and OUTP. 

The initial value of X, used in subroutine FCT and OUTP, specified 
by PRMT(l), Not necessarily the same as RO. (The observing 
station need not be assumed to be the surface of the earth.) 

This stands for "x and RO"; its initial and final values are 
entered as input data to specify the range of the independent vari- 
able X over which the ray path is to be calculated by ERFRAC. 

This is the current value of x minuts its initial value, which is its 
value at the observing station; in OUTP, 

Rectangular coordinates of the point on the reversed ray path as 
computed in the step immediately preceding the cimrent one; in 
OUTP. 


XWGB See RWGB. 

y This is the rectai^lar coordinate y in a coordinate system whose 

center is at the center of the earth and whose y-axis is parallel to 
a line drawn due east in the horizon plane of the station which ob- 
serves the ray being traced. See Fig. D-1. 

Y A four- element one dimensional array containing the elements 

listed below: 

Y(l) The variable name used for y in ERFRAC, 

Y(2) This is the first derivative of y with respect to x: dy/dx. 

It is also used in initialization procedures in ERFRAC to transfer 
the initial value of x, PRMT(l), from the main program to sub- 
routine FCT, where it is called XO. 

Y(3) The variable name used for z in ERFRAC. 

Y(4), This is the first derivative of z with respect to x; dz/dx. 
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YDEV, ZDEV 


YNOM, 

YOLD 

Y2SQ 

Y4SQ 

YWGB 

z 

ZDEV 

ZNOM 

ZOLD 

ZWGB 


The y (east) and z (north) components of the displacement of the 
current point on the ray path from a point at the same value of x 
on the projection of the reversed radio ray as observed at the 
observing station; in OUTP. 

ZNOM The y (east) and z(north) components of a point on the projection 

of the reversed radio ray as observed at the observing station, at 
the current value of x. 

See XOLD, 

Y(2) * Y(2). 

Y(4) * Y(4), 

See RWGB. 

This is the rectangular coordinate in a coordinate system whose 
center is at the center of the earth and whose z-axis is parallel to 
a line drawn due north in the horizon, plane of the station which 
observes the ray being traced. See Fig. D-1. 

See YDEV. 

See YNOM. 

See XOLD, 

See RWGB, 
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APPENDIX B 

THE DIFFERENTIAL EQUATIONS FOR THE RAY PATH 


Note: This appendix is reprinted from the Goddard Space Flight Center Report 
X- 520-69-30, pages 120 to 124, in which the rectangular coordinate axes were labeled 
with X eastward, y northward, and z along the observing station's vertical, z being the 
independent variable. But the IBM subroutines DRKGS and DHPCG use X as the inde- 
pendent variable, and I thought it wise to be consistent in this usage throughout the 
present program, ERFRAC. Consequently, as used in ERFRAC, the formulas reprinted 
in this appendix are modified by chaning z to x, x to y, and y to z. See also section 5 of 
the present paper which explains how the two second-order equations derived in Appendix 
B are replaced by a system of four first-order equations for use in ERFRAC, 

Those variables used in Appendix B, expressions B24 to B29, reprinted here, which 
are relabeled in ERFRAC, are as follows: 


Original 

Name 

Name 

ERFRAC 

Meaning 

DNOMZZ 

DENOM 

Eqn. B26 

MUX 

MUY 

B m/ 9 y 

MUY 

MUZ 

B /i/B z 

MUZ 

MUX 

B |ii/B X 

xz 

Y(2) 

B y/B X 

XZSQ 

Y2SQ 

(3y/3x)^ 

XZZ 

DERY(2) 

B^ y /fe x^ 

XZZNUM 

DY2NUM 

Eqn. B-25 

YZ 

Y(4) 

3 z/B X 

YZSQ 

Y4SQ 

(Bz/B x) ^ 

YZZ 

DERY(4) 

B^ z x^ 

YZZNUM 

DY4NUM 

Eqn, B-27 


THE DIFFERENTIAL EQUATIONS FOR THE RAY PATH 

Fermat's Principle states that the path taken between two points by a ray is that 
one such that, as compared with neighboring paths, the optical path length, the line 
integral of the index of refraction aloi^ the path, has a stationary value. The optical 
.path length, OPL, for any given path between points A and B, is given by the path integral 


OPL = r /ids 


(Bl) 


where is the index of refraction for the wave, in general a function of position in space, 
and ds is an element of length along the path. 
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In rectangular coordinates, 


M M (x, y, z) 


d s = (d + d + d z^)!/^ 


(B2) 


= (1 + x2 + y2)l/2 


d z - d 


where x = x(z) and y = y(z) are parametric equations of the ray path; 


and 


Xj = 3 X/3 z, y s 3 y/3 z; 


g = 1 + x2 + y2 


The optical path length can now be written: 


OPL = C fx (x, y, z) g ‘/2 d 2 _ r (x, y, z, 


X,. yj dz 


where 


F = /Li g*/2 


(B3) 


(B4) 

(B5) 


(B6) 


This integral has a stationary value, according to Fermat’s Principle, and it can be 
shown (see for example Matthews and Walker 1964, page 309) that this implies that the 
integrand function F satisfies the Euler Lagrange equations: 


iZ - ^ F 

3 x d z 3 x 


(B7) 


3F _ _d_ 3 F 
3 y d z 3 y^ 


(B8) 


The left side of Equation B7 is: 



The right side develops as follows: 


= 1/2 /J. (2 X ) = M 

B X z z 


(BIO) 


d/d z (BF/Bx^) = B/B z (B F/B xj + d x/d z B/3 x (3 F/B x^> 

+ d y/d z 3/B y (B F/B x^) + d x^/d z B/B x^ (B F/B x^) 

+ B y/B^ B/B y^ (3 F/B x^) = ‘1’ + ‘2’ + ‘3’ + ‘4’ + ‘5’ (Bll) 

where 

■1’ = g->/2 

‘2’ = x^ g-l/2] 

‘3’ = y^ l>ty 

‘4’ = x^^ [/X g*‘/2 + ^ (- 1/2) g‘3/2 (w x^)] (B12) 

= x^^ [m g"‘/2 - M x^ g'3/2]^ 

‘5’ = y„ [pLX^ (- 1/2) g-3/2 2 yj 

-3/2 

= x^y^y,, g 

When these results are substituted into Equation B7 and it is multiplied through by g'"'^ 
it becomes: 


g = M X + u 


+ 


y, + M x^ - M x^ 


- /2 X y y 


(B13) 


Rearranging terms: 


/X (1 - x2/g) - M x^ y^ g-i y^^ = g - M, x^ - M, x^ - x^ y^ 


^ /2-x g - (Mz + Xx + My y^) x^. 


(B14) 
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a = fjL 

b = M x2 g-l 
C = M Vj g“‘ 
d = M g“‘ 

e g 

f = My g ■ 

® + M, + Aty y^. 

Then Equations B14 and B15 become: 

(a - b) = d y^^ = e - h x^ 

and 

- d x^^ + (a - c) y^^ = f - h y^. 

The solution is: 

(e - h x^) - d 

(f - h y ) a - c 

(a-c)(e-hx ) + d(f-hy ) 

X = - Z ^ 

(a - b) - d (a - b) (a - c) - d2 

- d a - c 


(B16) 


(BIT) 

(B18) 

(B19) 
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(B20) 


a - b e - h X 

Z 

- d f - h I 

(a - b) (f - h y^) + d (e - h x^) 

(a - b) (a - c) - d^ 

For the special case where the ray path is entirely in the x plane, y(z) = 0, and from 
B16: 


c = 0 f = 0 
d = 0 h = /i.^ + x^ 

Then Equations B19 and B20 give: 


(B21) 


and 


a(e-hx)+0 e-hx 
2 ' _ 2 

"" “ (a - b) a - 0 ■ a - b 


S - Cm, + (g - x2-) - x^ 

M - At x2 g*l M (1 - x2/g) 


^ (a - b) (0 - O') + 0 
(a - b) a - 0 


:z 0 


For computer programming in the general case, we define: 

A r MU 

M = MU * XZSQ/G = MUOVG * XZSO 
C = MU * YZSQ/G = MUOVG * YZSQ 
D = MU * XZ » YZ/G = MUOVG * XZ * YZ 
E = MUX * G 
F = MUY * G 

H = MUZ + MUX * XZ + MUY * YZ 
MUOVG = MU/G; XZSQ = XZ * XZ; YZSQ = YZ * YZ 


(B22) 


(B23) 


(B24) 
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Then from Equations B19 and B20, we have; 


XZZNUM (A - C) * (E - H * XZ) + D * (F - H * YZ) (B25) 

DNOMZZ = (A - B) * (A - C) - D * D (B26) 

YZZNUM = (A - B) * (F - H * YZ) + D * (E - H * XZ) (B27) 

XZZ = XZZNUM/DNOMZZ (B28) 

YZZ = YXXNUM/DNOMZZ (B29) 
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APPENDIX C 


VALUES OF DYDXO AND DZDXO 

The initial values of y(2) = dy/dx and y(4) = dz/dx for an upward projected ray at 
the observing station can be written down in terms of the observed altitude, e and 
azimuth, 4> , of the ray, by inspection of Fig. C-1: 



d X = d s sin 0 


(Cl) 


d y = d s cos 9 sin 4> 


(C2) 


d z = d s cos 9 cos <j> 


(C3) 


d y/d X = cos 9 sin 0/sin 9 = cot i9 sin 0 


(C4) 


d z/d X = cos 9 cos 0/sin 9 = cot 9 cos 0 


(C5) 
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In ERFRAC terminology, 9 = ALTOBS, and <^ = AZOBS, and we have 


and 


where 


(2) = CXITALT * DSIN (AZOBS), 

(C6) 

(4) = COT ALT * DCOS (AZOBS), 

(C7) 

COTALT = DCOTAN (ALTOBS) 

(C8) 
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APPENDIX D 


COORDINATES OF RAY-WEDGE BASE INTERSECTION 


Here we derive formulas for the rectangular coordinates of the point at which a 
straight-line projection of the initial direction of the reverse ray at the observing station 
reaches a given altitude HWGB. This altitude is called simply h in this derivation. 

The actual reversed racjio ray will reach altitude HWGB at very nearly the same 
point, since its total bending never exceeds a very few minutes of arc, and one minute 
of arc at even 300 km distance corresponds to a sideways displacement of less than 
0.1 km. 

The problem is to determine the point, P, at which a straight line intersects a 
sphere of radius r^ + h, where tj, is the radius of the earth. (See Fig. D-1). Since the 



Figure D-1. Geometry for Calculating Coordinates of 
Ray-Wedge Intersection 


line goes through the point (x^, 0, 0), which is the location of the observing station, the 
equations for it can be written: 


(x - XplZ-t = y/m = z/n, 


(Dl) 


where (I, m, n) are the direction cosines of the line. The equation for the sphere is: 
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x2 + y2 + z2 = (r^ + h)2 = r2, 


(D2) 


where we define 


r 


+ h 


Expressions for y and z from Eq. Dl, when substituted into Eq, D2, give: 

+ (m (x - X|j)/'6')2 + (n (x - ■x^)/'l)'^ = 
or 

x2 (1 + mV^2 ^ r,2/^2^ - 2 X (m^ + nVl^) 

+ x2 (m2/^2 ^ ^ r2 

or 

x^ - 2 Xg (m^ + n^) X + (m^ + n2) x^ = r^, 

where we have used the fact that the direction cosines satisfy the relation 

= 1 . 

If, for brevity, we introduce 


Eq. D5 becomes: 


p = = 1 - 


x^ - 2 Xg p X + (p x^ - r^) = 0 

The more positive of the two roots is: 


X = (2 Xj, p + (4 x2 p2 - 4 p x2 + 4 r2)l/2/2 


= P + <^’‘0 " ’‘o P + r2)i/2. 


(Consideration of Fig. D-1 shows that we want this one of the two intersections 
line theoretically makes with the sphere.) 


(D3) 

(D4) 

(D5) 

(D6) 

(D7) 

(D8) 

(D9) 

(DIO) 
(Dll) 
which the 
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The substitution of Eqs. D3 and D8 for r and p gives, after simplification; 


X = Xg (1 - + -t ((r^ + h')^ _ (1 _ l'^') x^')!/^ 


(D12) 


As a simple check on the correctness of this equation, we try it for the case '6 = 1, which 
corresponds to a ray traveling vertically upward. The equation then gives the intersection 
point at X = Cq + h, which is correct. Also, for the case ^ = 0, the ray travels along a 
station horizontal, so that the answer given by the equation, x = i;, , is again correct. 

In ERFRAC terminology, Eq. D12 is: 


XWGB = (1. - LSQ) * X0 + L * DSQRT ((R0 + HWGB) * (R0 + HWGB-) 
- (1. - LSO) * X0* X0) 


(D13) 


where LSQ is L * L. 

The variables RO and HWGB are among those read in near the beginning of sub- 
routine FCT and are therefore available in storage. The direction cosine L is the cosine 
of the angle between the reverse ray and the vertical and is therefore the sine of the 
observed ray altitude angle: 


-t = sin 9 


(D14) 


A previous computer step gave to DERY(4) this value: 


DERY (4) = DSIN (ALTOBS), 


(D15) 


SO that it could be transmitted to subroutine OUTP. Therefore here we merely set: 


L = DERY (4) 


(D16) 


..and 


LSQ = L * L 


(D17) 


Once the value of x at the intersection is obtained by Eq. D13 the values for y and 
z follow from Eq. D1 : 


y = (m/'t) (x - Fg) 

(D18) 

z = (n/'t) (x - Fjj) 

(D19) 
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In terms of the observed altitude e and azimuth 0 , 


y = (cos 0 sin <p/sin 0) (x - 


= cot 0 sin <p (x - r^); 


(D20) 


z =: (cos 0 cos 4>/sin 0) (x - 

(D21) 

= cot 0 cos 0 (x - r^) 


The coefficients of (x - r ) in these two equations are, of course, just the values of the 
derivatives dy/dx and dz/dx for the reversed ray at the ground. These have previously 
been evaluated and stored as y(2) and y(4). (See App. C.) Therefore 


YWGB = Y (2) ♦ (XWGB - R0-) 


(D22) 


and 


ZWGB = Y (4) * (XWGB - R0-) 


(D23) 


(In ERFRAC, y is the eastward coordinate and z is the northward coordinate.) 
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APPENDIX E 


DIRECTION COSINES OF THE ELECTRON WEDGE GRADIENT 


The electron wedge base plane in ERFRAC is taken as tangent to an earth- centered 
sphere of radius (r^ + h^), at the point P where the projection of the reversed ray 
reaches altitude h^, as shown in Fig. D-1 of Appendix D. The "electron wedge gradient," 
as the term is used here, lies in the base plane, at an angle AZWG with the plane defined 
by the station vertical (the x axis) and the point P, as shown in Fig. E-1. 



Figure E-1. Geometry of Electron Wedge Gradient 


Consideration of this figure shows that the direction cosines of the wedge gradient 
with respect to the x, y, and z axes established at the ground station, are respectively 
given by 
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AWG = - cos AZWG sin 6 ^ 

BWG = sin AZWG cos AZOBS 

+ cos AZWG cos 0JJ sin AZOBS 

CWG = - sin AZWG sin AZOBS 

+ cos AZWG cos 6 ^ cos AZOBS 


Since the coordinates of the wedge base point P are XWGB, YWGB, and ZWGB, cos 
0Q in Fig, E-1 is given by: 


cos = COSTH = XWGB/RWGB, 


where 


RWGB = (XWGB2 + YWGB2 + ZWGB^)!/^ 


Also, then, sin is given by: 

sin = SINTHT = (1 - C0STHT2)J/2. 

Finally, we introduce 


and 


SINAZW = DSIN (^AZWG^) 


COSAZW = DC»S (^AZWG^) 


The formulas for the direction cosines then become: 


AWG 3 - CXJSAZW * SINTHT 
BWG = SINAZW * COSAZO 

+ COSAZW ♦ COSTHT ♦ SINAZO 

CWG 3 - SINAZW * SINAZO 

+ COSAZW * COSTHT * COSAZO 
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Unless the observed ray arrives at a lower angular altitude than is normal in Minitrack 
observations, the great circle angle between the station vertical and the wedge base 
point P {6 in Fig. E-1) will be small, perhaps two or three degrees. Then the geo- 
graphical "azimuth of the wedge gradient, , measiared easterly from north, will be 
approximately given by: 


4 >^ AZOBS + AZWG, (El) 

or 

AZWG - <^0 - AZOBS (E2) 


If AZOBS is specified, this equation enables one to determine the value of AZWG to use 
t simulate a wedge whose gradient is oriented at a given azimuth angle (approximately), 
from north. However, it may actually be of more immediate interest to specify the 
angle between the plane containing the projected reverse ray and the station zenith and 
the direction of the wedge gradient, and this is AZWG itself. 
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APPENDIX F 

FORMULAS FOR LOBS AND MOBS 


Inspection of Fig, C-1 shows that the direction cosines of the reversed ray, with 
respect to the y (east) and z (north) axes, respectively, are given by: 


LOBS = cos 0 sin <p 


(FI) 


MOBS = cos ^ cos 


(F2) 


The quantities Y(2) and Y(4) were previously calculated according to formulas C-4 
and C-5 of Appendix C, as follows: 


y (2) = cot 0 sin 4> 


(F3) 


y (4) = cot 0 cos <j> 


(F4) 


Since cos 0 = (sine) (cote), our formulas F-1 and F-2 can then be written: 


LOBS = cos e sin 4> - sin 0 cot 0 sin <p - y (2) sin 0 


MOBS = cos e cos 4> = sin 0 cot 0 cos <p = y (4) sin 0 


As part of the initialization procedures, the value of sin e is passed from subroutine 
FCT to subroutine OUTP as DERY(4). Therefore, we finally have, in ERFRAC variables: 


LOBS = dery (4) * y (2) 


MOBS = dery (4) * y (4) 
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APPENDIX G 



PERFO: 


rCE DATA FOR DRKGS AND DHPCG 


Data on CPU (central processing unit) time charged for the various trials of 
DRKGS and DHPCG are tabulated below. The table shows the number of runs that were 
made at each value of IHLF used. 




DRKGS 

DHPCG 

Jobs 

of 

69.218 

IHLF 

0 

2 

4 

6 

0 

1 

3 

5 

7 

No. of Runs 

2 

1 

1 

1 

1 

1 

1 

1 

1 

*Total CPU Time 

3.45 sec 

4.51 sec 


Jobs 

IHLF 

4 

5 

6 

7 

5 

6 

7 

of 

69.220 

No. of Runs 

6 

12 

6 

6 

12 

6 

6 


♦Total CPU Time 

57.35 sec 

37.82 sec 


*Note: CPU Time used for program compilation and link steps, about one second, is not included. 


In order to calculate the computer time required by each of the subroutines, 
DRKGS and DHPCG to make a given run, for purposes of comparison, I here assume 
that, usii^ a given subroutine, a run with IHLF increased by one takes twice as much 
time. 


This seems reasonable to me, since when the integration step size is cut in half, 
twice as many steps must be computed to cover the specified range of X (X = 0 to X = 7.) 
Let r and p be the time required for DRKGS and DHPCG, respectively, to make a single 
run with step size corresponding to IHLF = 5. Then using the assumption stated above, 
the data for the jobs of 69.218 give: 


or 


and 


or 


2 (r/32) + 1 (r/8) + 1 (r/2) + 1 (2 r) = 3.45, 


(G2-1) 


r = 1.28 seconds, 


(G2-2) 


1 (p/32) + 1 (p/16) + 1 (p/4) + 1 (p) + 1 (4 p) = 4.51, 


(G2-3) 


p = 0.84 seconds 
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The data for the jobs of 69.220 give: 


6 (t/2) + 12 r + 6 (2 r) + 6 (4 r) = 57.35, 


(G3-1) 


or 


r = 1.12 seconds , 


(G3-2) 


and 


12 (pj + (2 p) + 6 (4p) = 37.83, 


(G3-3) 


or 


p = 0.79 seconds 


While the data for the two jobs do not lead to exactly the same results for r or for 
p, they are in reasonable agreement, and indicate that the computer time required by 
DHPCG to complete a single run with a given step size comes out significantly less than 
that required by DRKGS, 34% less by the first calculation and 29% less by the second. 
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APPENDIX H 


ELECTRON DENSITY IN THE ELECTRON WEDGE 


A side view of an electron wedge region as used in ERFRAC, is included in Fig. 
D-1, near the top of the figure. An enlarged view is shown in Fig. H-1. 





/ 


TO CENTER / 

OF EARTH ^ 


Figure H-1. Electron Wedge Geometry 


The direction cosines of the wedge gradient with respect to the x, y, and z axes 
are called AWG, BWG, and CWG. (See Appendix E.) The overall thickness, TWG, of the 
wedge at an arbitrary point (x, y, z) is taken to be the thickness, TWGO, at the base 
point plus an additional amount which depends upon the distance SWG, that the point 
(x, y, z) is from the base point measured parallel to the "wecfee gradient." If Rg is the 
vector drawn from the base point to the point (x, y, z), and s is a unit vector in the direc- 
tion of the wedge gradient, then SWG is clearly given by the dot product of Rg and s. 

Since the components of s are AWG, BWG, and CWG, this can be written; 


SWG = AWG (X - XB) + BWG (Y - YB) + CWG (Z - ZB), 


(H2-1) 
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where XB, YB, and ZB are used here to represent the coordinates of the wec^e base 
point. The local thickness of the wedge is then given by: 


TWG = TWG0 + (SWG) (SLOPWG) 


(H2-2) 


Let HWG be the perpendicular distance of the arbitrary point above the base plane 
of the wedge, as shown in Fig, H-1. Now we assume that the electron density due to the 
wedge depends upon HWG in accordance with the following formula: 


newg = newmx ♦ exp [1 - (1 - . 


(H3-1) 


where q, defined by 


q = 2 (hwg/twg) - 1. , 


(H3-2) 


varies linearly with HWG from q = -1, at the base of the wedge, to q = +1, at its top. 

The function specified by equation (H3-1) symmetrical about q = 0; it is shown in 
Figure H-2 for q = 0 to q = +1, It closely resembles a parabola except that at each end 
of its range it has a point of inflection and a short but smooth tail by which it approaches 



Figure H-2. Assumed Electron Density Distribution 
in the Electron Wedge 
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zero tangentially. The author tried using a simple parabolic distribution, which ap- 
proaches zero at q = ±1 with a finite slope, but the Integrating subroutine was unable to 
get past this point of discontinuity in the gradient of total electron density. He came 
across the function used in equation (H3-1) while browsing in a book on elementary 
topology.) 

NEWGMX is equation (H3-1) is the maximmn electron density in the wedge, which 
is assumed to be constant, independent of wedge thickness. This assumption is defen- 
sible in ERFRAC because of the fact that the ray being traced will never pass through 
the region near the vertex of the wedge, regardless of what orientation parameters are 
specified for the wedge. (See Appendix E.) 

The parameter HWG can written in terms of known variables by noting that it 
is the dot product of the vector R g, introduced above, and a unit vector r which is per- 
pendicular to the base plane, and therefore parallel to a line drawn from the center of 
earth to the base point (XWGB, YWGB, ZWGB). If RWGB is the length of this line, 

HWG is then given by; 


HWG = (XWOB/RWCB) (X - XWGB) + (YWGB/RWGB) (Y - YWGB) 


+ (ZWGB/RWGB) (Z - ZWGB) 


(H6-1) 


We can now derive formulas for the space partial derivatives of the electron density. 
Consider first the partial derivative with respect to x. Using lower case abbreviations 
of the variable names appearing in the formulas above, and a prime (') to represent 
partial differentiation with respect to X, we have, from equation (H3-1) 


n = exp (1 - (1 _ q^)"') 


(H7-1) 


n' = n (1 -q2)-2 (_ 2 q) q' 


(H7-2) 


= - 2 n q (1 - q2)"2 q’ 


Now from equation (H7-2), 


q' = 2 (h'/h t'/t)/t 


= 2 (h' - h t'/t)/t 


(H7-3) 


From equation (H6-1) 


h' = Xg/R^, 


(H6-1) 
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from equation (H2-2), 


t' 


= s' p. 


(H6-2) 


and from equation (H2-1), 


(H6-3) 


Using these results in equation (H7-3) and (H7-2), we obtain: 

q' = 2 (x^/Rj^ -ha p/t)/t 


n' = - 2 n q Cl - q2) ^ (2) (x^/Rj^ -ha p/t)/t 


(H9-1) 


or 


B n/B X = + 4 * newg * (1 - q^)”^ (hwg * awg * slopwg/twg (H9-2) 

- xwgb/rwgb/rwgb)/twg 

Because of the symmetry in equations (H2-1) and (H5-1), the formulas for Bn/By 
and for dn/dz can immediately be written down by replacing XWGB and AWG with 
YWGB and BWG, to get Bn/By, or with ZWGB and CWG to get Bn/Bz , If we define: 


= 4 * newg *g(l. - q*q)*(l.-q*q)twg 

(H8-1) 

and 


B = HWG * SLOPWG/TWG, 

(H8-2) 

then the three derivatives can be written: 


B N^B X = (B * AWG - XWGB/RWGB), 

(H8-3) 

B Ny By = A^ (B * BWG - YWGB/RWGB-), 

(H8-4) 

and 


B NyB z = Ajj (B * CWG - ZWGB/RWGB). 

(H8-5) 


where we have used in place of n to represent electron density. 



(H5-2) 


n = n sin^ q 

m ^ 

" = "n, - ‘5^') 


Differ entiatii^ this: 


n' = - 2 q n q' 

n'= 2n^(sinqcosq)q' 

Now from equation H3-1, 

q' = 2 (h'/t - h t'/t2). 

= 2 (h' -h t'/t)/t 


From equation (H5-1), 


(H5-3) 


(H5-4) 


165 



IT" MyP-ENDIX I 


INDEX OF REFRACTION AND ITS DERIVATIVES 


Lawrence et al (1964) give as an approximate expression for the index of refraction, 
in the ionosphere: 


M = 1 -bN /o)2, 

^ e 


(II -1) 


where 


b = e^/ 2 m. 


( 11 - 2 ) 


Here, N^ is the number of free electrons per cubic meter, w = 277f is the angular 
frequency of the radio wave, c is the charge on an electron, is the permittivity of 
free space, and m is the mass of an electron. When o> = Inf and the values of the con- 
stants in MKS units are substituted, the formula Il-l becomes: 


/li = 1 - 40.3 N /f2 


(11-3) 


In ERFRAC, N, is in units of 10^“ electrons/m® and f is in units of 10® Hz, Us- 
ing these units, equation 11-3 becomes: 


M = 1 - 0.403 N^/f2 


( 12 - 1 ) 


For convenience in computation, the refractivity, N, is defined as follows: 


N = (/i - 1) X 10® 


(12-2) 


From equation (12-1), 


N = - 4.03 X 10* N /f2 

e 


(12-3) 


Since /j- = 1 + 10‘® N, it is clear that 


B m/B X = 10"S B N/B X. 


(12-4) 


From equation (12-3), 


B N/B X = - (4.03 X 10Vf2) B N^/B x 


(12-5) 
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Then equation (12-4) becomes: 


3 /n/B X = - (.403/f^) 3 x 
Now using equation (H8-3) of Appendix H, 

3 /^/3 X = - (.403/f^) \ * AWG - XWGB/RWGB), 

where Aj, is given by equation (H8-1). 


If now we define the new variable A as follows: 


then we have 


A = - (,403/f^') A„. 


MUXWG = 3 /i/3 X = A (B * AWG - XWGB/RWGB) 


MUYWG = A (B * BWG - YWGB/RWGB) 


MUZWG = A (B * CWG - ZWGB/RWGB) 


(13-1) 


(13-2) 


(13-3) 


(13-4) 


(13-5) 


(13-6) 


(13-7) 
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APPENDIX J 


INPUT DATA FORMAT 


As originally designed and written, eight cards of input data are required for 
ERFRAC, as follows: 


Card 

No, 

Format 

Variables, In Order 

(Refer to List of Variables in Appendix ,) 

1 

4F10.0 

Note: All distances are to be specified in meters. 
RO; XMRO.tart ; XMRO^nd ; INITIAL STEP SIZE. 

2 

D8.2 

PRMT(4). (Upper bound on truncation error. Enter in the form 

3 

F5.0 

bbl.D-mn, where bb represents two blanks and mn is the desired 
exponent, such as 09.) 

FREQ, (In megahertz.) 

4 

2F10.0 

ALTOBS; AZOBS. (Both in degrees.) 

5 

4F10.0 

SECSUN; HCHP; HEMAX; NEMAX 

6 

5F10.0 

NEWGMX; HWGB; TWGO; SLOPWG; AZWG. (Enter AZWG in 

7 

4F10.0 

degrees.) 

DERY(l); DERY(2); DERY(3); DERY(4). 

8 

FIO.O 

DXPRNT. 
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